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Abstract 


Research  in  multiscale  methods  has  recently  flourished  with  the  help  of  ever- 
improving  computer  technology.  These  developments  enable  computational 
physics  methods  to  challenge  many  of  the  fundamental  limitations  of  continuum 
mechanics  with  larger  atomistic  simulations  and  sophisticated  hybrid  atomistic- 
continuum  methods.  The  foundation  of  most  hybrid  methods  presently  lies  in 
the  judicious  application  of  kinematic  constraints  between  regions  of  atoms  and 
regions  of  continuum  finite  elements.  This  juxtaposes  atomic  and  continuum 
force  fields  and  introduces  an  interface  along  which  atoms  and  nodes  are 
unnaturally  constrained.  The  constraint  is  necessary  to  establish  compatibility  of 
displacements  across  the  interface.  This  report  is  divided  into  three  sections. 
Section  1  reports  on  an  investigation  of  finding  sources  of  numerical  error  due  to 
this  unphysical  constraint.  Bounding  estimates  on  the  numerical  error  are 
derived  using  summation  rules  of  a  classical  interatomic  potential  and  the 
geometry  and  periodicity  of  the  molecular  structure.  A  previously  observed 
inverse  relation  between  element  size  and  interface  error  is  demonstrated,  and 
additional  numerical  experiments  are  presented.  In  section  2,  a  literature  review 
of  a  technique  that  can  potentially  eliminate  this  error  is  presented.  The  review 
covers  efforts  in  engineering  for  composite  materials  rooted  in  a  firm 
mathematical  basis  for  the  so-called  asymptotic  expansion  homogenization 
method  (AEH).  The  homogenization  method  is  used  as  a  framework  for 
developing  a  multiscale  system  of  equations  in  elasticity,  also  in  section  2.  In 
section  3,  AEH  is  used  as  a  framework  for  developing  analytical  multiscale 
formulations  for  frozen  atoms  at  the  small  scale  and  continuum  mechanics  at  the 
large  scale.  The  Tersoff-Brermer  type  II  potential  (Brenner,  D.  W. 

Physical  Review  B.  Vol.  42,  no.  15,  pp.  9458-9471, 1990;  Tersoff,  J.  Physical  Review 
Letters.  Vol.  61,  no.  25,  pp.  2879-2882, 1988)  governs  the  atom  interactions,  and 
hyperelasticity  governs  the  continuum.  A  quasi-static  assumption  is  used 
together  with  the  Cauchy-Born  approximation  to  enforce  the  gross  deformation 
of  the  continuum  on  the  positions  of  the  atoms.  This  makes  the  atomistic 
equations  linear.  The  two-scale  homogenization  method  establishes  coupled 
self-consistent  variational  equations  in  which  the  information  at  the  atomistic 
scale,  formulated  in  terms  of  the  Lagrangian  stiffness  tensor,  feeds  the  material 
information  to  the  continuum  equations.  Analytical  results  in  one  dimension  are 
shown. 
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1.  Estimating  Numerical  Error  in 
Atomistic-Continuum  Computational 
Methods  in  Graphene 


1.1  Overview 


The  underlying  premise  of  the  finite  element  method  (FEM),  when  applied  to  solving  prob¬ 
lems  in  computational  mechanics,  is  that  the  material  in  question  is  a  continuous  medium 
infinitely  divisible  into  smaller  continuous  components.  Developments  in  the  basic  sciences 
such  as  chemistry  and  physics,  however,  have  disputed  this  assumption  long  before  the  no¬ 
tion  of  FE  was  even  first  conceived.  In  most  applications  at  macroscopic  length  scales,  it  is 
useful  to  think  of  materials  as  a  system  of  continuously  distributed  mass.  Yet  today,  with  the 
development  of  advanced  devices  such  as  nano- /microelectromechanical  systems  (N/MEMS) 
whose  features  and  characteristics  challenge  the  most  fundamental  assumptions  of  contin¬ 
uum  mechanics,  it  is  clear  that  traditional  FEM  cannot  be  applied  without  additional  con¬ 
siderations.  Moreover,  the  numbers  of  atoms  in  these  devices  warrant  the  determination  of 
properties  that  require  computations  on  a  scale  unreachable  by  atomistic  methods  alone. 

Methodologies  for  linking  a  continuum  to  an  atomistic  domain  can  be  found  in  the  litera¬ 
ture  as  early  1971  [1]  in  which  the  treatment  of  the  continuum  is  predominantly  atomistic 
in  nature.  Finite  element  methods  were  later  employed  by  Mullins  and  Dokainish  [2]  using 
a  numerically  decoupled  domain  approach  with  spatially  overlapping  atomistic  and  contin¬ 
uum  regions  in  which  the  information  from  each  region  is  fed  into  the  other  via  boundary 
conditions.  A  review  of  these  methods  can  be  found  in  Cleri  et  al.  [3].  Among  these  early 
analytic  and  computational  studies,  frequent  issues  regarding  the  treatment  of  the  interface 
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arose  which  were  primarily  handled  through  creative  use  of  kinematic  constraints. 

More  recently  Tadmor  et  al.  [4]  develop  an  entirely  FE-based  formulation,  the  so-called 
quasicontinuum  method.  Local  and  non-local  formulations  are  used  to  discriminate  between 
atoms  in  regions  of  low  interest  and  high  interest.  The  transition  between  the  two  regions  is 
handled  again  through  kinematic  constraints  and  the  mesh  is  adaptively  refined  in  regions 
of  large  atom  motion. 

The  so-called  handshaking  or  coupling  of  length  scales  (CLS)  method  [5]  adds  an  additional 
level  of  sophistication  with  a  third  region  modeled  using  the  tight-binding  (TB)  method 
which  captures  information  about  the  electronic  degrees  of  freedom.  There  are  three  regions 
of  interest:  (1)  an  FE  region  modeled  using  elasticity  theory,  (2)  a  molecular  dynamics  (MD) 
region  using  a  classical  potential,  and  (3)  the  TB  region.  The  coupling  at  interfaces  sepa¬ 
rating  these  regions  is  accomplished  through  kinematic  constraints  at  the  FE/MD  interface 
and  chemical  constraints  at  the  MD/TB  interface. 

A  generalized  scaling  approach  is  developed  in  coarse-grained  molecular  dynamics  (CGMD) 
[6]  to  better  handle  the  propagation  of  waves  through  the  atomistic-FE  interface  and  the  FE 
far  field.  A  coarse  graining  procedure  is  used  to  “merge”  the  atomistic  degrees  of  freedom 
from  Hamilton’s  equations  of  motion  to  the  smaller  number  of  FE  nodal  degrees  of  freedom. 

The  common  theme  in  the  previously  mentioned  investigations  is  to  directly  connect  a  con¬ 
tinuum  region  to  an  atomistic  region  either  through  kinematic  or  statistical  constraints.  The 
benefit  of  these  types  of  approaches  is  the  removal  of  complexity  in  regions  of  the  domain 
where  detailed  atomic  resolution  is  unnecessary  by  replacing  a  large  number  of  atomic  de¬ 
grees  of  freedom  with  a  smaller  number  of  element/nodal  degrees  of  freedom.  This  type  of 
approach,  however,  does  not  ensure  compatibility  at  the  interface  and,  as  a  result,  so-called 
“ghost  forces”  may  occur  [7]. 
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Friesecke  and  James  [8],  as  part  of  an  investigation  to  derive  a  scheme  that  passes  an  atomistic 
energy  to  a  continuum  energy,  showed  that  the  sum  of  the  energies  of  in  the  elements  is 
greater  than  (less  negative)  the  total  energy  of  the  original  body  before  discretization.  In 
fact,  the  difference  between  the  sum  and  total  energies  is  proportional  to  0(1/7*).  This 
provides  an  indirect  explanation  for  the  ghost  forces.  The  point  of  departure  in  this  work  is 
in  the  explicit  calculation  of  the  energy  error  based  on  geometry  arguments  of  the  molecular 
structure  of  the  material. 

There  are  two  primary  ways  of  classifying  the  forces:  those  due  to  the  initial  discretization 
or  modeling  error  and  those  due  to  subsequent  deformations  of  the  mesh  that  stem  from  the 
discretization.  The  aim  of  this  section  of  the  report  is  to  quantify  the  first  type  of  error  in 
terms  of  energy  and  to  specify  their  origins  from  a  classical  potential.  We  will  only  briefly 
discuss  the  second. 

In  practice,  creative  constraints  have  been  applied  near  the  interface  region  using  transition 
schemes  to  reduce  the  error  [2,7].  However,  a  systematic  means  of  approximating  this  error 
is  still  unavailable.  The  objective  of  this  section  is,  therefore,  to  make  a  first  step  in  this 
direction  by  developing  an  approach  where  this  error  can  be  estimated  directly  from  the 
classical  potential  and  the  structure  of  the  molecular  geometry.  This  has  application  to 
problems  where  the  mechanical  state  of  the  system  depends  on  the  total  energy.  Examples 
may  include  problems  of  phase  transition  and  local  material  instability.  The  present  potential 
is  specifically  suited  for  semiconductor  lattices  with  a  two-atom  basis,  namely  carbon,  and 
accounts  for  two-  and  three-atom  effects  [9]. 

In  section  1.2.  some  preliminary  definitions  are  given;  in  section  1.3,  a  discussion  of  the 
sources  of  numerical  error  in  atomistic-continuum  computational  schemes  is  presented.  In 
section  1.4,  a  description  of  the  actual  classical  potential  used  in  the  present  calculations 
is  described  followed  by  estimates  of  the  energy  error  and  discussions  of  the  results.  Some 
additional  considerations  for  more  general  studies  will  then  be  presented  in  section  1.5, 
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followed  by  final  remarks  in  section  1.6. 


1.2  Problem  Definition 


Consider  a  lattice  of  periodically  spaced  atoms  in  the  n-dimensional  space  ft  E  Rra  with 
na  atoms  and  nb  bonds  connecting  the  nearest  neighbor  atoms.  The  atoms  are  initially 
undistorted  and  at  zero  temperature.  The  plane  of  atoms  can  be  fully  described  by  a  set  of 
primitive  cell  vectors  (e*)  and  additional  translation  vectors  (p j)  if  there  is  more  than  one 
atom  in  each  basis  cell.  For  carbon  in  the  form  of  a  single  graphite  sheet,  graphene,  there 
are  two  atoms  in  the  basis  (j  =  1)  and  therefore  only  one  translation  vector. 

The  sheet  of  atoms  is  permitted  to  deform.  The  deformation  is  measured  by  the  deformation 
gradient  F.  Let  the  deformation  u  :  ft  -»■  with  gradient  F  =  have  an  energy  density 
given  in  the  general  form 

-i  ft b 

w=-y>(F),  (i) 

b 

where  b  runs  over  all  the  bonds  (nb)  in  ft,  and  the  ^  coefficient  gives  the  energy  density  per 
atom.  Assume  that  $  is  defined  for  all  F  that  are  real  2x2  matrices  with  det  F  >  0  and 
that  $  exhibits  the  standard  features  of  lattice  invariant  deformations  in  ft. 

The  energy  density  expression  can  take  any  general  form  of  a  cluster  potential.  Therefore, 
the  summation  can  involve  bonds,  bond  angles,  and  higher  bond-order  effects.  A  specific 
case  study  is  made  in  section  1.4  which  requires  a  more  exact  description  for  the  potential. 
However,  the  discussions  here  are  generally  applicable  for  any  classical  potential. 

For  the  sake  of  discussion,  assume  that  the  energy  in  equation  (1)  can  be  decomposed  into 
the  sum  of  an  energy  due  to  a  deformation  and  a  non-zero  reference  equilibrium  energy.  The 
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(2) 


reference  equilibrium  energy  is  based  on  the  initial  positions  of  the  atom  nuclei,  X,  i.e., 

-j  rib  1  nb 

W  =  -£>(F)  +  -5>£(x), 

na  b  na  b 

where 

*,(F  =  I)  =  0.  (3) 


Through  a  change  of  variables,  equation  (2)  can  be  rewritten  as  a  function  of  the  atom 
displacements, 

1  _n&_  1  Ub_ 

w=-y>(p(u))+-y>?(x).  (4) 

na  "  na  ^ 

The  assumption  that  the  energy  is  separable  in  equation  (2)  is  nontrivial.  In  general  classical 
potentials,  the  ansatz  is  highly  nonlinear  and  behaves  poorly  in  the  traditional  partial  dif¬ 
ferential  equation  sense.  It  possesses  at  least  one  singularity,  is  nonconvex  and  noncoercive. 
The  simplification  introduced  in  equation  (2)  will  be  justified  next  through  a  bilinearization 
of  the  problem  that  removes  this  difficulty  at  the  expense  of  generality. 


In  light  of  the  small  deformation  assumption,  we  can  take  a  bilinearization  (harmonic  ap¬ 
proximation)  about  the  equilibrium  lattice  configuration.  This  involves  taking  the  Taylor 
series  expansion  of  the  energy  as  follows, 


W  =  W|F=I  + 


dW 


dF 


F=I 


:  (F  —  I)  + 


(5) 


F=I 


When  the  system  is  in  equilibrium,  the  first  derivative  is  zero.  Omitting  higher  order  terms, 
this  leaves 


By  definition, 


w=  W|P=I  +  i  (F  -  I)r  : 


d2W 

WdF 


(F  -  I). 


(6) 


F  -  I  =  Vu. 


(7) 
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Substituting  equation  (7)  into  equation  (6)  gives 


w  =  W|F=I  +  -  (Vu) 


1  ,T 


dFdF 


:  (Vu) . 


F=I 


(8) 


Notice  that  the  reference  equilibrium  energy  occurs  naturally  in  the  expansion.  This  justifies 
the  earlier  assumption  of  decomposability. 


As  a  consequence  of  the  small  deformation  assumption  and  the  resulting  harmonicity  in  W, 
we  can  assume  u  €  H1.  Or  more  conventionally,  one  can  say  u  €  <S  for  a  given  set  of 
conditions  of  u  on  fiu,  where 


S  =  {u|u  €  H\  u|nu  =  u}.  (9) 

Furthermore,  define  V  C  S  such  that 

V  =  {w|w  e  H\  w|nw  =  0}.  (10) 


Then,  for  a  given  set  of  prescribed  displacements  u  on  the  boundary  Qu,  the  stable  configu¬ 
ration  of  the  atoms  is  one  which  minimizes  the  total  potential  energy 


n(u)  =  inf 


f  WdQ  —  f  f  •  vdCl  —  f 
i/  q  v  J  nQ 


t  •  vdT 


(11) 


where  v  G  S  are  the  trial  functions,  f  is  the  body  force  per  unit  volume  and  t  is  the  prescribed 
surface  tractions  per  unit  area  on  Clg.  Substituting  equation  (8)  into  equation  (11)  gives, 

dW 

'  lF= 


Il(u)  —  fl0  =  inf 


t  •  vdT 


(12) 


Associated  with  equation  (12)  is  a  bilinear  symmetric  form  a(-,  •)  and  inner  products  (-,  •) 


and  (•,  -)r  where 


a(w,  u)  = 


d2W 

dFdF 


:  (Vu)dO, 


(f ,  u)  =  f  f  •  udft, 
Jn 

(t,  u)r  =  [  t  •  udr. 

Jng 


(13) 

(14) 

(15) 
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Then,  the  functional  in  equation  (12)  has  the  form, 

'1 


Il(u)  -  II0  =  inf 


tt(v,v)  -  (f,v)  -  (t,  v)r 


(16) 


with  the  equivalent  weak  form  of  equation  (16)  given  by:  Find  u  €  S  for  every  w  €  V  such 
that 


o(w,u)  =  (w ,  f)  +  (w,t)r- 


(17) 


The  finite  element  analogue  of  (17)  is  based  on  the  conventional  discretization  of  displace¬ 
ments  and  corresponding  gradients, 

Tie 

=  YlWiUh  (1§) 

i 

Vu"  =  ]Tb  iUi,  (19) 

i 

where  i  =  1,  ...,ne  are  the  nodes  in  the  mesh,  Wi  are  the  weighting  functions,  and  Bj  are  the 
gradients  of  the  weighting  functions.  The  finite  dimensional  approximations  of  S  and  V  are 
denoted  by  Sh  and  Vh,  where  Sh  C  S  and  Vh  C  V. 


1.3  Discretization  Error 


Classical  potentials,  such  as  empirical  interatomic  potentials,  are  phenomenological  in  nature 
and  do  not  explicitly  model  the  effects  of  the  electron  density.  However,  tests  have  shown  that 
in  studying  strain  energies,  particularly  in  the  context  of  mechanical  deformation,  classical 
potentials  give  equally  reliable  results  compared  to  ab  initio  calculations  (Figure  1)  for 
homogeneous  deformation  under  the  present  assumption  of  small  strains  (<10%).  Despite 
the  myriad  sources  of  error  that  can  stem  ostensibly  from  the  use  of  a  classical  potential, 
the  aim  of  this  work  is  to  study  the  error  that  arises  due  to  the  introduction  of  the  finite 
element  discretization. 
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Figure  1 .  Deformation  energy  density  comparisons  for  homogeneous  stretch¬ 
ing  of  graphene  in  the  plane.  Smooth  line  denotes  the  classical 
potential  and  dots  denote  computed  ab  initio  values. 


The  internal  energy  must  now  be  expressed  in  an  approximate  finite  element  form.  Discretize 
Q.  with  a  set  of  elements  Sleh  such  that 

^  =  (2°) 

e 

which  has  the  associated  minimization  problem  given  by 

n'V )  -  n*  =  inf  [£afc(v\v*)  -  ( f,Vy  -  (t,  ,  (21) 

V*  |_z 

where 

a‘(w\“‘)  =  ^  (Vu'")dnft,  (22) 


(f.uY  = 

f  f  •  u hdQh, 

(23) 

Jsih 

II 

f  tu hdTh. 

(24) 

J  n* 
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As  in  equation  (17),  the  associated  weak  form  is  given  as:  Find  €  Sh  for  every  €  Vh 
such  that 


ah(wh ,  uft)  =  ( wh ,  f)h  +  (wfe,  t)p. 


(25) 


It  remains  to  more  carefully  define  the  bilinear  form  ah(-,  •)•  From  equation  (4),  the  internal 
energy  is  written,  as  a  consequence  of  the  discretization,  as 

1  ne  nb  ne  nb 

w‘  =  -£X>(F(u))  +  ;r£5>‘(x>>  (26) 

Tia  •  ,  'I'd  u 

ib  ib 

where  h  runs  over  all  the  elements  and  b  runs  over  all  the  bonds  within  the  element,  n\. 
Then,  define  €  Sh  as  the  solution  that  satisfies  equation  (25)  for  Wh  given  by 

1  Tib  -t  ne  Ub 

vv*  =  ^-£^(F(“h))  +  -  E£^<x)-  (27> 

Ua  b  a  i  b 

subject  to  the  prescribed  loads  and  displacements  on  the  boundary. 


The  driving  motivation  of  atomistic-continuum  problems  is  the  reduction  of  the  number  of 
degrees  of  freedom.  However,  through  the  introduction  of  the  finite  element  approximation, 
two  sources  of  ghost  forces  arise.  Notice  that  these  can  now  be  elucidated  from  equations 
(4),  (16),  (21),  and  (27).  The  first  of  these  will  be  described  briefly  here.  The  second  will 
be  expounded  in  section  1.5.  The  first  occurs  at  initial  equilibrium,  i.e.,  u  =  0,  where  the 
difference  in  total  potential  energies  between  exact  and  approximate  is  the  difference  in  the 
reference  lattice  energies.  Let  the  lattice  of  all  bonds  in  Q,  be  denoted  Ub  and  those  in  Q.h  be 
denoted  L\.  Then,  the  energy  difference  is  given  by 


n-nA  = 


n0-nh0 


na 


Tle  na 


£X>«<X> 


£  «x) 


(28) 
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where  b  runs  over  all  bonds  that  cut  across  element  boundaries,  i.e.,  all  bonds  that  connect 
atoms  in  two  different  elements.  The  discretization  of  leads  to  a  set  of  disjoint  sets  of  atoms 
that  causes  the  two  terms  in  equation  (28)  to  be  unequal  generally.  The  only  instance  when 
they  are  equivalent  is  when  ne  =  1  and  na  =  n®,  or  when  Cb\C £  G  0  and  no  discretization  is 
performed. 

Equation  (28)  is  an  equilibrium  energy  form  of  the  ghost  force  error  and  is  due  to  a  more  pre¬ 
dominant  assumption  that  stems  from  a  discretization  error  as  evidenced  by  the  summation 
being  over  all  bonds  excluded  in  the  meshing.  Bond-counting  arguments  and  summation 
rules  can  then  be  used  to  characterize  and  quantify  this  error  near  the  interface. 

There  are  two  observations  from  the  energy  difference  in  equation  (28).  First,  it  is  nonzero 
at  equilibrium.  Approximating  Q  with  flh  artificially  introduces  internal  surfaces  into  the 
problem,  effectively  truncating  the  communication  among  the  atoms  (Figure  2). 

Although  this  effect  grows  asymptotically  smaller  for  very  large  elements  (relative  to  atom 
spacings),  the  error  is  more  pronounced  as  the  element  size  becomes  comparable  to  atom 
spacings.  This  situation  occurs  frequently  in  atomistic-continuum  computational  methods 
due  to  the  kinematic  constraints  needed  to  enforce  local  compatibility  between  the  atoms 
and  the  elements.  Second,  at  small  scales  equation  (28)  is  neither  translation  nor  rotation 
invariant,  a  characteristic  of  isotropic  continua,  because  of  mesh-dependence  at  small  scales. 
Figure  3  illustrates  this  by  showing  that  the  energy  in  element  A  is  unequal  to  the  energy  in 
element  B  despite  both  elements  having  the  same  perimeter  and  areas.  This  is  the  nature  of 
discretizing  a  problem  that  contains  an  already  discrete  set  of  atoms. 
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(a)  Original  graphene  sheet  with 
bonds  denoted  by  straight  lines  and 
C-atoms  denoted  by  black  dots. 


(b)  Fragment  of  an  unstructured  fi¬ 
nite  element  mesh  used  to  approximate 
atomic  degrees  of  freedom  with  nodes 
(open  circles). 


Figure  2.  Element  boundaries  act  as  truncation  lines  that  cut  off  the  com¬ 
munication  among  atoms  in  different  elements.  The  result  is  the 
artificial  increase  in  the  total  energy. 


1.4  Example:  Graphene 


We  now  attempt  to  quantify  the  error  for  the  specific  example  of  graphene  using  the 
Stillinger- Weber  classical  potential  given  by  Stillinger  and  Weber  [9]  and  Abraham  and 
Batra  [10].  It  has  been  documented  to  give  reliable  lattice  binding  energy  and  atom  spacing 
for  graphite.  The  energy  is  given  by 

T(iibrw)x{r(ij))x(nik)),  (29) 

i  j  i  j  k 

i<j  i^j  ^3rk 

where  (p  and  'tp  are  the  respective  two-  and  three-atom  terms,  x  is  the  short-range  cut-off 
function,  ( i ,  j ,  k )  are  atom  indices,  and  the  density  is  retrieved  by  dividing  by  the  number 
of  atoms  W  =  W/na.  The  two-  and  three-atom  potentials  and  cut-off  function  are  given 
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Figure  3.  The  total  energies  of  elements  A  and  B  are  not  equivalent. 


respectively  by 


A(Br(ij)  ~  Kr(«)  -  a)  *] 

0 


^,(r(*i)’r(  **)) 

x(Hij))  = 


eX 


r(b) 

-r(y) 


-  cos(0*) 

r(ik) 


< 


exp  [7 (rW)  -  o)  *] 
0 


r  <  a 
r  >  a 


5 


r  <  a 
r  >  a 


(30) 

(31) 

(32) 


where  A,  a,  A,  #*,  7  are  material-specific  constants,  and  the  symbol  r^)  denotes  the  vector 
originating  from  atom  i  and  terminating  at  atom  j.  This  material  was  chosen  specifically 
because  of  its  covalency  as  a  semiconductor  and  the  absence  of  significant  long-range  forces. 
For  graphene,  only  short-range,  nearest-neighbor  forces  need  to  be  considered. 


We  introduce  a  local  coordinate  system  centered  on  the  graphene  primitive  cell.  Any  straight 
line  that  passes  through  the  cell  can  be  exactly  defined  by  the  set  ( ee,9e ).  Figure  4  depicts 
a  straight  line  cutting  through  the  primitive  cell  and  the  associated  local  coordinate  system. 
The  straight  line  is  an  idealization  of  an  element  boundary  cutting  through  the  atomistic 
region.  It  represents  the  line  across  which  atoms  are  not  permitted  to  “communicate.”  As 
a  heuristic  upper  bound  on  the  magnitude  of  the  error,  (ee,6e)  can  be  permuted  so  that  the 
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Figure  4.  Local  coordinate  system  and  cutting  line  that  represents  fictitious 
element  boundary. 


line  cuts  through  the  largest  number  of  bonds  per  unit  length.  This  corresponds  to  the  set 
(ee,  6e)  =  (qe2, 0)  for  |  <  q  <  §.  The  symmetry  of  the  hexagonal  rings  also  implies  there  are 
six  rotations  that  give  equivalent  results. 


The  interesting  result  here  is  the  upper  bound.  Each  bond  has  an  energy  equal  to  the  pair 
potential  value  of  6  for  the  undeformed  nearest  neighbor  atom  spacing  of  r0.  For  the  (ee,  9e ) 
upper  bound  parameters  previously  given  on  a  periodic  cell,  two  bonds  are  cut  per  r0\/3 
distance.  Therefore,  the  energy  error  density  per  unit  area  of  the  graphene  sheet  is 

n  -  nA  =  c\,  C  <  (33) 

h'  ~  rj 3 

Through  numerical  experiments,  the  efficacy  of  this  upper  bound  can  be  tested  for  increasing 
element  sizes.  Figure  5  shows  several  example  calculations  for  varying  ee  and  6e.  The 
data  sets  marked  0i,62,03,6A.6o  respectively  denote  (ee,0e)  =  (4e2,0),  (|e2,0),  (|e2,^), 
(2e2’ik)>  (|e2,§),  anc*  (3e2-0)-  The  scatter  at  smaller  element  sizes  shows  the  breakdown 
of  rotational  and  translational  invariance  near  the  atomistic  region  as  previously  discussed. 
The  scatter  is  also  attributed  in  part  to  the  imprecise  measure  of  bond  density  when  only 
a  few  bonds  are  involved  at  very  small  element  sizes.  Notice  that  the  magnitude  of  the 
error  for  small  h  is  as  high  as  3  eV/A  with  a  mean  value  of  approximately  2.2  eV/A.  The 
equilibrium  interatomic  distance  in  graphene  is  approximately  1.42  A.  The  standard  lattice 
binding  energy  is  —7.44  eY  per  carbon  atom.  Then,  in  the  extreme  case  where  there  is 
poor  mesh  selection,  the  energy  error  can  exceed  4  eV/atom.  For  larger  element  sizes,  the 
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Figure  5.  Numerical  tests  on  the  discretization  error  (|Wft|  =  C). 

coefficient  of  the  error  approaches  a  continuum  limit  where  invariance  sets  in,  and  the  error 
of  the  total  energy  density  is  even  smaller  due  to  the  inverse  dependence  on  h. 

A  lower  bound  can  be  generated  likewise  by  choosing  the  best  approximating  scenario  of  the 
fictitious  element  boundary.  If  we  assume  that  any  line  that  cuts  directly  through  an  atom 
affects  neither  the  atom  nor  the  bonds  emanating  from  it,  then  the  lower  bound  is  always 
zero  since  a  line  can  always  be  oriented  in  a  periodic  graphite  lattice  to  cut  only  through  a 
line  of  atoms.  An  example  of  this  is  the  set  ( ee,0e )  =  (962,0)  for  9  =  |  or  |.  Specifically, 
the  coefficient  is  found  to  be 

C  >  (34) 

3  Tq 
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This  calculation  suggests  that  refinement  is  not  always  desirable  to  atomic  scales  when 
employing  finite  element  methods,  an  observation  that  is  counterintuitive  to  conventional 
h-adaptive  methods.  A  means  of  systematically  correcting  these  effects  is  still  currently 
unavailable. 


1.5  Nontrivial  Deformation 


The  second  source  of  ghost  forces  comes  from  the  first  terms  in  equations  (4)  and  (27),  i.e., 


1  Tla  -i  Tla 

a(u,u)  -  aft(u\uft)  =  —  y>.(F(u))-— £>„(F(u*)) 

Tla  7la 


(35) 


r,„  ,  d2w 

X(VU)  :  dFdF  IM 

_  f  (VuA)  : 

J»  '  dFdF 


:  (Vu )dn 
:  (Vu h)dnh. 


(36) 


F=I 


Despite  the  common  expansions  about  F  =  I  in  both  integrals  in  equation  (36),  the  La- 
grangian  stiffness  terms  are  still  unequal  because  of  the  differences  in  the  implied  summa¬ 
tions  over  the  requisite  atoms  (and/or  bonds).  More  precisely,  the  number  of  bonds  that  are 
counted  in  W  is  unequal  to  the  number  in  Wh.  Thus,  standard  numerical  analysis  techniques 
cannot  be  applied  directly  without  significant  modification. 


In  the  limit  as  h  — »  oc  and  r0  — >  0,  it  has  been  shown  that  this  atomistic  error  exhibits  a 
boundedness  proportional  again  to  O  (£)  [8].  Furthermore,  under  these  assumptions,  one 
can  show  that  the  stiffness  terms  are  equivalent  in  the  limit,  and  the  convergence  rate  of 
equation  (36),  the  square  of  the  energy  norm,  is  bounded  by  some  function  of  the  element 
size  h.  Or  in  standard  notation 

||e||n<c/ifc-n+1||u||fc+1,  (37) 

where  e  =  uh  —  u,  c  is  some  independent  constant,  k  is  the  complete  order  of  the  piecewise 
smooth  polynomial  used  in  the  interpolation  and  n  is  the  order  of  the  derivative  in  the 
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energy  expression,  i.e.,  the  appropriate  metric  space.  For  continua,  these  types  of  rigorous 
estimates  are  well  established  [11].  However,  for  problems  involving  resolvable  atomic  details 
relative  to  finite  element  size,  the  precise  magnitudes  of  error  have  not  yet  been  ascertained 
in  a  generalized  way  as  in  continua.  Such  developments  are  potentially  useful  for  error 
estimator/corrector  schemes  for  implementation  into  existing  computational  methods. 

Under  limited  distortions,  however,  simple  estimates  as  in  section  1.4  can  again  be  obtained 
for  deformed  configurations.  We  presently  consider  the  case  where  the  in-plane  shearing  de¬ 
formation  of  the  atoms  coheres  to  the  Cauchy-Born  rule  [12].  Furthermore,  the  deformations 
are  limited  to  the  class  of  all  uniform  shears  in  which  the  coordination  number  of  the  atoms 
stays  fixed.  That  is,  the  interatomic  spacings  do  not  change  so  that  new  bonds  do  not  form 
and  existing  bonds  do  not  break.  The  shear  is  characterized  by  the  deformation  gradient 

(  1  a  0  ^ 

F  =  0  10  ,  (38) 

V°  0  1). 

where,  based  on  the  cut-off  function  in  equations  (29)  and  (32),  and  all  relevant  material 
parameters  for  carbon,  we  choose  0  <  a  <  0.0617.  The  deformed  configuration  of  graphene 
is  depicted  in  Figure  6.  Under  homogeneous  deformation,  the  three  vectors  a,  b,  and  c  are 
the  only  vectors  needed  to  describe  the  bond  orientations  and  lengths  in  the  entire  sheet  by 
appropriately  tiling  the  plane.  The  Cartesian  components  of  the  deformed  bonds  are  given 

by 


The  restriction  on  a  makes  this  a  first  order  approximation.  Pair  interactions  dominate 
over  the  higher  order  triple  terms  in  equation  (29).  Using  the  earlier  arguments  for  the 
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Figure  6.  Homogeneous  shear  of  graphene. 

undeformed  configuration,  the  UB  and  LB  estimates  for  the  error  coefficients  are 

=  —  [</>(o)  H-  +  2^(a,  b)  +  2ip(a,  c)  +  2,0(b,  c)] , 
ru 

CLB  =  —  [2V»(a,c)  +  2^(a,b)], 

Li, 

where 

rj/  =  —  6\/3a;  +  9a2,  =  \/3r0, 

z 

and  o  ■  |a|  and  6  =  |b|. 

The  change  in  the  error  estimates  are  shown  vs.  the  extent  of  deformation  in  Figure  7.  The 
variations  are  in  the  meV  range.  Further  considerations  must  be  undertaken  to  evaluate 
the  error  under  larger  deformations  where  changes  in  coordination  and  higher  order  atom 
interactions  occur. 


(40) 


(41) 


1.6  Closing  Remarks 


The  area  of  computational  multiscale  method  development  remains  a  fertile  area  of  research. 
It  attempts  to  answer  many  of  the  most  difficult  questions  of  science  at  the  nanoscale  by 
leveraging  the  understanding  of  continua.  In  this  report,  we  have  proposed  a  first  order 


Figure  7.  UB  and  LB  estimates  for  the  line  density  of  the  error  as  a  function 
of  deformation  (|V\4|  =  C). 

characterization  of  the  energy  error  in  undeformed  and  limited  deformed  configurations  of 
graphene  using  classical  energy  expressions.  In  the  present  carbon  system  where  the  force 
field  is  dominated  by  pairwise  interactions  of  atoms,  the  findings  show  that  the  error  in¬ 
troduced  by  discretization  may  lead  to  significant  errors  in  the  energy  under  certain  mesh 
refinements. 
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2.  A  Literature  Review  of  the  Asymptotic 
Expansion  Homogenization  Method 


The  first  part  of  section  2  briefly  describes  several  key  developments  in  homogenization. 
Their  deficiencies  in  being  unable  to  provide  local  information  motivates  the  subsequent  dis¬ 
cussion  of  the  so-called  asymptotic  expansion  homogenization  (AEH)  approach.  In  section 
2.2,  various  types  of  related  and  available  homogenization  methods  are  described  that  pos¬ 
sess  the  homogenization/localization  capability  to  handle  complex  linear/nonlinear  behavior 
encountered  in  continuum  and  phenomenological  engineering  problems.  Section  2.3  intro¬ 
duces  the  mathematical  literature  from  which  AEH  finds  its  basis.  Other  names  found  in 
the  literature  for  this  approach  include  asymptotic  homogenization,  mathematical  homoge¬ 
nization,  and  classical  homogenization.  All  are  based  on  the  fundamental  premise  that  the 
displacement,  velocity,  or  temperature  field  is  representable  by  an  asymptotic  series.  To 
avoid  confusion,  the  method  will  be  referred  to  as  AEH,  whereas  classical  homogenization 
techniques  will  refer  to  homogenization  methods  based  on  classical  mechanics  principles.  In 
section  2.4,  application-oriented  efforts  for  AEH  for  linear  problems  are  summarized.  Sec¬ 
tions  2.5  and  2.6  describe  inelastic  and  nonlinear  problems  in  the  AEH  literature.  The 
method  is  completely  demonstrated  for  a  linear  elastic  example  in  section  2.7.  And  final 
conclusions  are  drawn  from  the  literature  review  in  section  2.8. 

2.1  Classical  Homogenization  Methods  and  the  Intro¬ 
duction  of  AEH 


Early  efforts  in  finding  the  effective  properties  of  composite  materials  frequently  approached 
the  problem  from  a  mathematical  point  of  view  to  predict  bounds  on  properties  [13-16].  The 
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simple  geometries  of  spherical  or  cylindrical  inclusions  permitted  the  use  of  simplified  ana¬ 
lytical  methods.  Though  these  investigations  reveal  much  about  the  mechanics  of  composite 
materials,  assumptions  on  the  geometry  of  the  microstructure  and  ordering  of  constituent 
properties  pose  limitations  on  their  applicability  to  advanced  high  performance  composites 
with  complex  microstructures. 

Recent  composite  technologies,  such  as  sophisticated  woven  fabric  composites,  employ  com¬ 
plex  patterns  in  the  weaving  of  reinforcement  fibers.  Thus,  bounds  or  other  approximations 
of  the  effective  properties  that  were  originally  designed  for  simple  microstructures  such  as 
regular  arrays  of  uniform  cylinders  or  uniformly  distributed  spheres,  are  inappropriate.  The 
same  difficulties  exist  in  predicting  other  constitutive  properties  such  as  conductivity,  per¬ 
meability,  diffusivity,  thermal  expansion,  and  the  like.  Furthermore,  though  a  repeating  cell 
may  still  be  identifiable  and  assumed  small,  the  assemblage  of  cells  within  the  true  global 
body  may  not  exhibit  homogeneity  in  the  necessary  statistical  sense  [17]  making  the  direct 
application  of  analytical  techniques  difficult. 

To  incorporate  the  details  of  the  weaving  patterns  of  the  bundled  fibers  and  the  depen¬ 
dence  of  microstructural  parameters  on  the  macro  properties,  some  researchers  have  used 
laminated  plate  analytical  continuum  methods  [18-23].  By  varying  the  boundary  conditions 
applied  according  to  displacements  or  tractions  consistent  with  the  theorems  of  minimum 
potential  and  complementary  energies,  bounds  for  the  effective  plate  stiffness  and  compliance 
were  obtained.  Unfortunately,  these  methods  are  either  inherently  one  dimensional  (1-D)  or 
two  dimensional  (2-D)  or  are  restricted  to  simple  microstructures.  They  employ  the  so-called 
mosaic  or  crimp  approximations  for  the  woven  fabric  geometry  that  allows  them  to  construct 
analytical  expressions  for  the  fiber  shape  and,  consequently,  the  constitutive  equations.  As 
such,  they  do  not  provide  insight  into  the  three-dimensional  (3-D)  physics  associated  with 
advanced  composites,  particularly  when  the  fibers  are  not  in  the  shape  of  simple  trigono¬ 
metric  functions  or  when  complex  boundary  interactions  dominate  the  internal  behavior  of 
the  material. 
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Other  efforts  appear  in  the  form  of  mechanics-of-materials  approaches  within  the  framework 
of  numerical  methods  for  determining  the  homogenized  properties  [24-28].  These  approaches 
are  useful  for  complex  3-D  microstructures.  They  typically  require  some  external  boundary 
conditions  to  be  applied  to  the  unit  cell  from  which  an  understanding  of  the  response  of 
the  structure  is  obtained  by  comparisons  with  Hooke’s  Law  or  energy  balance  principles. 
Averaging  operations  of  the  phenomenological  models  over  the  unit  cell,  referred  to  as  ho¬ 
mogenization.  yield  the  homogenized  properties.  Details  of  the  unit  cell  are  smeared  away 
in  favor  of  a  single  set  of  properties  that  characterizes  the  inhomogeneous  material  in  an 
average  sense.  The  mechanics-of-materials  approaches  perform  homogenization  for  arbitrar¬ 
ily  shaped  microstructural  features.  The  AEH  approach,  on  the  other  hand,  can  handle 
arbitrary  shapes  but  also  enables  one  to  return  to  the  microlevel  details  in  the  unit  cell  (i.e., 
localization)  based  on  the  global  solution.  Whereas  the  mechanics-of-materials  approaches 
provide  homogenized  properties,  the  AEH  approach  provides  homogenized  properties  as  well 
as  localized  information. 

As  composite  manufacturing  techniques  grow  in  sophistication,  so  too  does  the  demand  for 
computational  and  analysis  methods  for  predicting  and  modeling  the  manufacturing  process. 
A  step  in  this  direction  is  computational  approaches  that  handle  the  important  multiscale 
issues  typically  associated  with  heterogeneous  materials.  Among  the  approaches  discussed 
in  this  study,  the  AEH  method  appears  viable  and  useful  for  determining  the  effective  prop¬ 
erties  of  composites  and  employing  the  computed  properties  in  subsequent  macrolevel  anal¬ 
yses.  The  approach  is  predicated  on  two  fundamental  assumptions:  that  the  displacements 
(or  other  primary  variable  such  as  temperature  or  flow  velocity)  can  be  characterized  by  an 
asymptotic  series  in  s  and  that  the  salient  features  of  the  microstructure  are  contained  in  a 
unit  cell  representative  of  the  periodicity  at  the  local  scale.  The  first  two  terms  of  the  series 
are  assumed  to  be  representative  of  the  respective  macro  and  micro  responses  (see  Figure 
8).  These  definitive  macro  and  micro  terms,  then,  subject  to  the  periodicity  constraints, 
provide  tangible  variables  from  which  to  extract  local  effective  micromechanical  properties, 
local  gradients,  and  the  overall  global  response.  The  approach  can  also  be  extended  to  ther- 
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(b)  Mathematical  illustration 

Figure  8.  The  body  in  e-space  is  the  realistic  representation  of  the  heteroge¬ 
neous  structure.  The  AEH  approach  isolates  the  micro  from  the 
macro  by  approximating  e  with  an  homogenized  body  in  X  and  a 
scaled  representative  unit  cell  in  Y. 

mal,  flow  or  structural  continuum  problems.  Its  large  range  of  computational  applications 
indicates  its  generality  for  a  broad  class  of  boundary  value  problems  (BVP). 

The  AEH  approach  has  also  been  studied  in  the  context  of  problems  involving  three  or 
more  scales  [29]  by  using  a  reiterative  scheme.  By  applying  two-scale  AEH  successively 
over  multiple  length  scales,  the  AEH  approach  can  be  extended  to  problems  involving  three 
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or  more  length  scales  with  little  modification  to  existing  formulations.  The  approach  has 
been  studied  mathematically  for  linear  heat  conduction  by  Bensoussan  et  al.  [29]  and  in 
engineering  analysis  for  linear  elasticity  in  Chang  and  Kikuchi  [30]. 

The  multiple  scale  abilities  notwithstanding,  the  AEH  approach  can  also  be  used  as  an 
approach  only  to  estimate  the  homogenized  properties.  By  deriving  partial  differential  equa¬ 
tions  (PDEs)  that  govern  the  influence  of  inhomogeneities,  no  restrictions  are  placed  on  the 
size  or  complexity  of  shapes  of  the  microstructure.  The  only  limitation,  therefore,  is  the 
significant  preprocessing  modeling  effort  needed  to  create  the  complex  geometries  associated 
with  modern  composites.  Thus  the  same  approach  may  be  applied  in  a  generalized  manner 
to  woven  fabric  composites,  knit-fiber  composites,  metal  matrix  composites,  and  the  like. 

Researchers  have  also  found  AEH  more  apt  for  extension  to  the  inelastic  regime  due  to  its 
ability  to  estimate  microlevel  information  and  continually  update  the  homogenized  macrolevel 
properties.  With  the  localization  capability  inherently  derivable  through  the  formulations, 
microlevel  information  such  as  yield  criteria  or  other  nonlinear  effects  and  state  variables 
can  be  updated. 

With  the  estimated  information  at  the  local  level,  global  effects  can  then  be  estimated  by 
re-homogenization  based  on  the  new  local  information.  The  new  global  information,  then, 
influences  the  global  homogenized  problem  in  a  circuitous  manner  [31-34].  Other  popular 
homogenization  approaches  [13-28]  do  not  permit  the  estimation  of  micro-  and  macrolevel 
information  simultaneously.  This  key  distinction  is  the  redeeming  quality  of  AEH  that  sets 
it  apart  from  other  homogenization  approaches. 

Conditions  encountered  in  the  analysis  of  composite  materials  are  usually  over  a  range  of 
temperatures  or  structural  behavior  making  them  nonlinear.  Traditional  micromechanical 
techniques  are  suited  best  for  conditions  resembling  linear  or  mildly  nonlinear  problems 
employing  many  simplifying  assumptions  on  the  geometric  shapes  or  constitutive  models 
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to  make  them  tractable.  Examples  of  restrictive  assumptions  may  involve  a  simple  mi¬ 
crostructure  geometry  [13-17]  or  linear  material  models  [18-22,24-26].  Their  ability  to  treat 
the  micromechanical  details  of  the  material  makes  them  useful  for  many  applications.  But 
the  majority  of  the  methods,  without  simplifying  assumptions,  are  normally  not  applicable 
to  nonlinear  situations.  Though  other  developments  such  as  those  in  nonlinear  continuum 
theories  have  been  extensive,  they  are  again  inherently  incapable  of  handling  the  complex  mi¬ 
cromechanical  shapes  and  the  aforementioned  multiple-scale  issues  associated  with  advanced 
composites. 

Though  investigations  in  the  literature  for  multiscale  problems  favor  the  AEH  approach  due 
to  its  fundamental  developments  in  functional  analysis,  other  homogenization/localization 
approaches  are  available.  However,  they  suffer  from  several  limitations.  The  approaches  and 
the  reasons  for  the  choice  of  AEH  are  described  in  the  next  section. 


2.2  Other  Homogenization/Localization  Methods 


This  study  employs  AEH,  whose  fundamental  approximation  is  that  the  primary  dependent 
variable  can  be  composed  by  the  superposition  of  a  smooth  global  solution  and  a  rapidly 
oscillating  local  solution.  The  sum  of  the  two  solutions  is  represented  by  the  first  two  terms 
of  an  asymptotic  series  in  e.  Based  solely  on  this  approximation,  the  underlying  goal  is  to 
develop  hierarchical  continuum  field  equations  for  the  coupled  multiscale  problem.  The  AEH 
method  is  not  the  only  approach  available  which  can  handle  the  complex  coupling  in  length 
scales  for  linear/nonlinear  applications.  Homogenization  approaches  for  linear  problems 
notwithstanding,  methods  of  homogenizing  material  properties  and  subsequently  localizing 
global  solutions  also  appear  elsewhere  in  the  literature.  The  present  state  of  progress  in 
these  areas  appears  to  be  in  the  early  mathematical  development  stages  or  without  sufficient 
generality  that  renders  them  useful  only  for  the  specific  applications  for  which  they  are 
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derived.  Few  of  these  methods  have  yet  to  be  used  routinely  in  engineering  applications  for 
practical  situations. 

A  detailed  discussion  and  survey  of  four  notable  approaches:  (1)  the  Fourier  Series  ap¬ 
proach,  (2)  the  Green’s  Function  approach,  (3)  the  Self  Consistent  approach,  and  (4)  the 
Subvolume  approach,  are  presented  by  Walker  et  al.  [35]  with  emphasis  in  the  analytical 
formulations.  Despite  the  seemingly  veritable  merit  in  the  approaches  and  the  subsequent 
references  therein,  few  researchers  appear  to  employ  these  methods  in  computational  me¬ 
chanics  elsewhere  in  the  literature.  This  is  in  contrast  to  the  sustained  growth  of  AEH 
methods.  The  ease  with  which  variational  equations  or  PDEs  can  be  formulated  via  AEH, 
unlike  the  methods  described  in  Walker  et  al.  [35],  explains  the  popularity  of  the  asymp¬ 
totic  approaches.  Though  the  present  study  emphasizes  mainly  AEH,  it  is  prudent  in  future 
investigations  to  scrutinize  the  methods  of  Walker  et  al.  [35]  closely  in  the  framework  of 
computational  mechanics. 

Moulinec  and  Suquet  [36]  employ  Fast  Fourier  Transforms  (FFT)  to  “pixelize”  complex 
geometric  microstructures  and  employ  a  superposition  of  displacements  to  develop  an  ho¬ 
mogenization  approach  to  study  elastic  and  inelastic  behavior  of  composites  via  FEM.  The 
FFT  circumvents  meshing  difficulties.  The  approach  for  nonlinear  situations  involves  a  step- 
by-step  integration  in  time  and  incorporates  the  exact  Green’s  function  of  the  linear  elastic 
and  homogeneous  analog  of  the  material.  Thus,  the  solution  is  expected  to  diverge  from  the 
true  behavior  of  highly  nonlinear  materials  where  the  Green’s  functions  depart  significantly 
from  elastic  and  homogeneous  conditions.  In  this  way,  the  approach  is  restricted  mainly  to 
the  linear  regime.  However,  it  shares  the  generality  of  the  AEH  approach  in  its  ability  to 
treat  a  large  variety  of  microstructures. 

Hou  and  Wu  [37]  integrate  FEM  with  an  AEH  approach  to  formulate  a  method  for  homog¬ 
enizing  arbitrary  heterogeneous  structures  not  limited  to  periodic  media  but  with  rapidly 
oscillating  microstructure.  The  developments  presented  are  for  elliptic  problems  which  in- 
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elude  elasticity  and  flows  in  porous  media.  The  fundamental  difference  between  this  effort 
and  other  efforts  in  AEH  is  in  the  description  of  the  base  functions.  The  approach  is  compu¬ 
tationally  expensive  due  to  the  solution  of  simultaneous  equations  required  to  obtain  these 
complicated  base  functions.  The  multiscale  base  functions  are  adaptive  to  the  local  prop¬ 
erties  to  account  for  the  refined  scale  of  the  heterogeneities.  Thus,  for  general  large  scale 
problems,  this  approach  becomes  computationally  prohibitive. 

Woo  and  Whitcomb  [38]  present  a  global/local  FEM  to  estimate  local  stresses  in  specific 
discretized  regions  using  a  refined  mesh  with  loads  kinematically  matched  at  the  subcell 
boundaries.  The  method,  however,  is  limited  primarily  to  linear  elastic  problems  due  to  its 
use  of  mechanics  of  materials  to  estimate  the  effective  properties  and  the  localized  behavior. 
In  contrast,  the  AEH  approach  does  not  employ  kinematic  conditions  when  establishing  a 
causal  link  between  the  micromechanical  behavior  and  the  global  response.  Instead,  the 
link  is  established  at  the  level  of  the  primary  variable,  e.g.,  the  deformation  for  structural 
problems,  via  an  asymptotic  series  approximation.  Thus,  unlike  the  AEH  approach,  the 
convergence  of  the  solution,  the  uniqueness  of  the  microlevel  stresses,  and  mathematical 
consistency  of  the  formulations  in  the  approach  of  Woo  and  Whitcomb  [38]  are  not  ensured. 

In  the  next  section,  the  background  and  related  efforts  in  AEH  are  described.  The  mathe¬ 
matical  basis  for  the  approaches  employed  in  this  study  can  be  found  in  a  large  breadth  of 
mathematics  literature.  These  are  described  next. 


2.3  Background  Literature:  Mathematical  Basis 


Although  asymptotic  methods  have  existed  for  many  decades,  the  application  of  an  asymp¬ 
totic  expansion  approach  for  heterogeneous  materials  is  relatively  new.  Some  early  efforts 
at  applying  perturbation  techniques  to  field  problems  in  a  strict  mathematical  context  are 
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scattered  through  the  literature  beginning  in  the  1960s  and  1970s.  The  literature  observed 
a  sharp  increase  of  published  efforts  in  the  late  1970s  and  early  1980s  in  AEH  as  engineers 
began  realizing  the  viability  of  the  method  for  transport  problems  and  as  the  maturity  of 
the  methods  influenced  others  in  the  mathematical  community.  The  increase  also  coincided 
with  the  release  of  two  authoritative  texts.  The  earlier  of  the  two,  by  Bensoussan  et  al.  [29], 
provides  specific  mathematical  definitions,  theorems,  and  principles  for  a  large  class  of  gen¬ 
eral  elliptic  systems.  In  addition  to  the  AEH  approach,  they  present  other  techniques  for 
the  treatment  of  heterogeneous  periodic  structures  in  general. 

The  second  text  by  Sanchez-Palencia  [39],  much  like  Bensoussan  et  al.  [29],  is  detailed  in 
its  mathematical  analysis.  Although  much  care  is  taken  in  presenting  the  equations,  very 
rarely  are  the  equations  actually  solved.  Yet  it  is  an  important  development  because  it 
provides  a  new  perspective  of  AEH  in  its  applicability  to  engineering  problems.  It  effectively 
demonstrates  the  utility  of  AEH  for  more  sophisticated  engineering  problems  by  deriving  the 
needed  equations  for  numerous  engineering  examples.  However,  both  Bensoussan  et  al.  [29] 
and  Sanchez-Palencia  [39]  favor  the  functional  analysis  aspects  of  AEH  and  provide  little 
discussion  on  how  to  solve  the  equations,  particularly  for  practical  problems. 

Other  AEH-related  efforts  in  the  literature  since  then  consistently  and  categorically  cite 
Bensoussan  et  al.  [29]  and  Sanchez-Palencia  [39]  as  their  mathematical  bases  for  subsequent 
novel  developments.  Other  mathematical  details  can  be  found  in,  for  instance,  Lions  [40], 
Duvaut  [41],  Oleinik  [42],  Tartar  [43],  Jikov  et  al.  [44],  Cioranescu  and  Paulin  [45],  and 
Bakhvalov  and  Panasenko  [46]. 

In  the  next  section,  the  literature  associated  with  the  AEH  approach  for  general  linear  compu¬ 
tational  mechanics  is  described.  Linear  computational  mechanics  encompass  low  Reynold’s 
number  flow,  heat  conduction/diffusion,  and  linear  elasticity. 
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2.4  Linear  Engineering  Problems 


The  majority  of  efforts,  as  documented  in  the  mathematics  literature,  focus  on  linear  prob¬ 
lems  primarily  due  to  the  ease  with  which  linear  problems  can  be  studied  via  functional 
analysis.  Despite  the  significant  number  of  publications  in  this  area,  discussions  herein 
are  limited  to  those  efforts  that  employ  AEH  in  an  engineering  context.  Engineering  is¬ 
sues  related  to  AEH,  namely  computational  issues  and  multiscale  equation  development,  are 
nontrivial  and  important  facets  to  understand  the  behavior  of  periodic  heterogeneous  media. 

Most  developments  are  presented  in  elasticity  due  to  the  myriad  range  of  engineering  situa¬ 
tions  that  can  be  covered  based  on  this  fundamental  elliptic  PDE.  Examples  include  elastic 
damage  and  fracture  mechanics.  Multidisciplinary  problems  are  also  considered  in  conjunc¬ 
tion  with  the  elasticity  problem,  namely  thermoelasticity,  which  requires  the  homogenization 
of  the  thermomechanical  properties.  The  homogenization  of  thermal  properties  can  also  be 
found  in  the  literature  for  efforts  in  solving  the  energy  equation  in  non-isothermal  fluid  flow. 

The  literature  here  is  arranged  according  to  areas  related  to  linear  transport  (fluid  and  heat 
flow)  or  linear  elasticity. 


2.4.1  Transport  Problems 

The  early  treatment  of  the  AEH  approach  can  be  found  in  areas  related  to  flow/transport  and 
heat  conduction  in  porous  media.  The  statistically  homogeneous  nature  of  porous  materials 
is  known  to  satisfy  the  periodicity  condition  required  in  the  AEH  approach.  A  straight¬ 
forward  presentation  of  AEH  for  flow  through  porous  media  is  shown  in  Keller  [47].  The 
homogenization  approach  is  employed  to  rederive  Darcy  s  Law  yielding  sets  of  equations  for 
the  two  length  scales:  micro  and  macro.  In  so  doing,  the  approach  avoids  common  heuris- 
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tic  methods  of  homogenization,  such  as  simple  volume  averaging,  by  introducing  discrete 
micromechanical  variables  via  higher  order  perturbation  terms.  The  velocity,  density,  pres¬ 
sure,  and  external  body  force  are  expanded  in  asymptotic  series  where  only  the  first  two 
terms,  representative  of  the  respective  macro  and  micro  variables,  are  carried  throughout 
the  formulations.  These  discrete  terms  introduce  traceable  quantities  and  provide  a  means 
of  ensuring  existence  and  uniqueness  of  solutions  to  the  multiple  scale  linear  problems. 

Citing  the  deficiencies  of  Darcy’s  Law,  namely  the  inability  to  handle  velocity  gradients  and 
boundary  layers,  Ene  [48]  presents  the  AEH  approach  for  the  Brinkman  equation.  Among 
other  more  complex  engineering  situations  also  considered  include  flow  in  a  fractured  porous 
medium  and  fluid-solid  interactions. 

Chang  and  Kikuchi  [49]  applied  the  homogenization  method  to  analyze  the  non-isothermal 
mold  filling  process  used  in  resin  transfer  molding  and  structural  reaction  injection  molding. 
The  approach  employs  a  doubly  porous  woven  fiber  preform,  the  first  scale  at  the  scale  of 
individual  filaments,  the  second  at  the  scale  of  bundled  fiber  tows,  and  the  third  at  the  scale 
of  the  global/macro  structure.  Stokes  flow  is  assumed  for  the  microscale  transport  where 
the  resistance  of  the  flow  provides  an  estimate  of  the  permeability  tensor  at  the  next  length 
scale.  The  tensor  is  then  used  in  Darcy’s  Law  which  can  then  be  homogenized  once  again  to 
provide  a  global  permeability  tensor  for  the  Darcy’s  Law  flow  assumed  to  govern  the  macro 
mold  filling  problem.  Hence,  the  homogenization  procedure  is  reiterated  to  model  the  pore 
structures  in  the  two  hierarchical  regimes. 

Flow  through  porous  media  is  primarily  a  linear  problem  because  the  Reynold’s  number 
is  typically  small.  Xo  AEH  investigations  appear  to  consider  nonlinear  inertial  effects  or 
convective  terms. 

Few  research  efforts  seem  available  for  AEH  solely  for  simple  heat  transfer  problems  apart 
from  the  mathematics  literature.  The  focus  here  is  on  engineering  problems.  The  exist- 
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ing  literature  in  this  area  appears  classifiable  into  two  specific  categories:  those  treating 
multidisciplinary  problems  and  those  considering  the  complex  microstructure  geometries  of 
certain  types  of  composites. 

Different  from  related  efforts  in  other  mechanics  fields,  thermal  properties  depend  little  on 
temperature  gradients  and  mostly  on  the  temperatures.  That  is,  it  is  more  common  to 
associate  thermal  conductivity,  k^,  as  a  function  of  temperature 

kij  =  ki,(T).  (42) 

Localization  of  the  primary  variable  is  typically  uninteresting  from  a  practical  point  of  view 
because  the  variable  is  assumed  to  adhere  to  an  asymptotic  series  whose  terms  after  the 
zeroth  order  are  very  small.  The  higher  order  terms  of  the  series,  representative  of  the 
oscillatory  temperature  variations  arising  from  the  periodic  distribution  of  inhomogeneities, 
are  small  compared  to  the  zeroth  order  behavior  of  the  temperature.  Therefore,  the  primary 
variable  at  the  local  level  will  appear  nearly  constant.  The  gradients  and  heat  flux  (or  strains 
and  stresses  in  solid  mechanics),  however,  provide  greater  information  and  show  significant 
nonuniformity  at  the  local  level.  In  light  of  these  observations,  it  is  expected  that  the  AEH 
heat  transfer  literature  is  limited  primarily  to  homogenization  with  few  considerations  of 
the  localized  information  because  the  local  temperatures  appears  constant  and  uninteresting 
due  to  the  dominance  of  the  first  term  of  the  expansion.  An  investigation  exploring  these 
observations  were  shown  by  Chung  et  al.  [50]. 

Thermal  problems  were  treated  first  analytically  by  the  early  mathematical  investigations. 
Many  of  the  relevant  references  are  shown  in  Bensoussan  et  al.  [29]  and  Sanchez-Palencia  [39]. 
In  contrast  to  the  governing  equations  in  other  mechanics  areas,  the  thermal  problem  is 
tractable  analytically  because  temperature  is  the  only  unknown  and  is  a  single  degree  of 
freedom. 

Aside  from  the  general  mathematical  considerations,  researchers  have  attempted  to  under- 
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stand  the  thermal  behavior  of  complex  composite  materials.  Woven  fabrics  are  employed 
in  the  manufacture  of  printed  wiring  board  substrates.  Dasgupta  and  Agarwal  [51]  use  the 
two-scale  AEH  approach  to  study  the  orthotropic  thermal  conductivity  of  plain-weave  fab¬ 
ric  composite  laminates.  An  analytical  subproblem  is  studied  to  augment  the  numerical 
approach.  Closed- form  expressions  are  proposed  to  estimate  the  effective  conductivity  of 
the  microscale  unit  cell.  The  analytical  results  are  compared  with  the  Finite  Element  (FE) 
solution  of  the  microscale  equations  and  limited  experimental  results.  The  analytical  aug¬ 
mentation  to  AEH  appears  useful  in  problems  where  the  solution  of  the  microscale  BVP 
must  be  avoided.  However,  two  concerns  are  evident  from  their  paper.  First,  the  solution  of 
the  microscale  BVP  provides  the  necessary  information  to  perform  localization  later  in  an 
integrated  mathematically  consistent  multiscale  analysis  whereas  their  analytical  approach 
provides  only  a  means  of  homogenization,  not  localization.  Second,  the  analytical  augmen¬ 
tation  is  tailored  for  specific  microstructure  geometries,  thus  removing  the  generality  of  the 
AEH  numerical  approach.  At  the  expense  of  the  ability  to  treat  general  microstructures, 
their  analytical  augmentation  to  AEH  removes  the  need  to  solve  the  microscale  equations. 
This  augmentation  approach,  which  does  not  appear  to  present  significant  computational 
cost  savings,  is  dubious  because  the  solution  of  the  microscale  equations  in  AEH  is  an  inte¬ 
gral  feature  of  making  the  approach  general  in  the  present  study. 

The  paper  by  Chang  and  Kikuchi  [49]  discussed  earlier  also  solves  the  energy  equation  in 
a  non-isothermal  mold  filling  problem  for  composites  process  modeling.  The  transient  heat 
conduction  problem  is  considered  assuming  a  linear  unchanging  conductivity  tensor  and 
surface  convection.  The  linearity  of  the  problem  leaves  the  heat  transfer  problem  uncoupled 
from  the  fluid  flow  portion  of  the  analysis.  It  is  of  the  opinion  of  other  investigators  that  the 
problem  must  be  considered  in  the  fully  coupled  context  for  accurate  modeling  [52]. 

Of  the  developments  to  date,  AEH  employed  in  the  study  of  elastic  materials  is  most  numer¬ 
ous  in  the  literature.  This  is  due  to  the  wide  range  of  problems  which  can  be  analyzed  by 
studying  linear  elasticity.  In  the  next  section,  a  brief  survey  of  the  literature  encompassing 
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elasticity  problems  for  AEH  is  presented. 


2.4.2  Elasticity  Problems 

The  bulk  of  the  developments  in  AEH  appears  for  the  elastic  problem.  This  is  due  to  the 
larger  number  of  applications  feasible  modeled  by  linear  elasticity  and  the  general  robustness 
of  formulations  for  continuum  mechanics  in  nonlinear  structural  mechanics  and  its  adapt¬ 
ability  to  multiscale  situations.  The  difficulty  of  problems  in  the  flow  and  heat  transfer  areas 
increases  significantly  even  for  small  departures  from  linearity.  Such  limitations  appear  not 
to  be  present  in  structural  mechanics. 

In  this  section,  the  literature  for  linear  elastic  problems  is  reviewed.  Emphasis  is  on  efforts 
that  employ  linear  elasticity  as  the  basis  for  applications  to  more  complicated  engineering 
problems  within  the  AEH  framework. 

Although  the  early  formulations  for  simple  elasticity  are  shown  precisely  in  Bensoussan  et 
al.  [29]  and  Sanchez-Palencia  [39],  interesting  and  novel  developments  arise  as  the  level  of 
complexity  increases  in  the  elastic  engineering  problem.  Lene  [53]  proposes  an  approach  for 
studying  damage  in  elastic  UD  composites  with  arbitrary  cross  section  where  the  damage 
origin  is  localized  to  fiber-fiber  and  fiber-matrix  interfaces.  The  damage  is  simulated  by 
a  displacement  discontinuity  applied  by  numerically  inserting  a  singular  element  with  zero 
thickness  (slipping  only)  in  the  damaged  region.  The  magnitude  of  the  discontinuity  is 
specified  by  a  Coulomb  friction  principle.  Throughout  the  analysis,  the  material  is  assumed 
elastic  and  additional  damage  evolution  laws  from  thermodynamic  theory  are  employed  to 
trace  the  damage  through  the  transient  problem. 

An  optimal  shape  design  methodology  is  introduced  by  Bendsde  and  Kikuchi  [54]  that  does 
not  limit  the  problem  to  equivalent  initial  and  final  design  topologies  and  avoids  FE  re- 
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meshing.  The  AEH  approach  is  employed  to  calculate  effective  elastic  properties  repeatedly 
to  satisfy  specified  design  requirements.  Upon  computing  the  optimal  distribution  of  material 
in  space,  an  anisotropic  mechanism  results.  The  mechanism  is  constructed  by  an  infimum 
of  periodically  distributed  small  holes  in  an  homogeneous  isotropic  material. 

Adaptive  mesh  refinement  to  improve  numerical  and  spatial  accuracy  is  considered  by  Guedes 
and  Kikuchi  [55]  for  effective  linear  elastic  composite  properties.  Convergence  and  error 
estimate  studies  are  also  undertaken  in  the  context  of  FEM.  The  AEH  approach  is  employed 
and  explicit  details  and  physical  interpretation  of  the  characteristic  function  or  correctors 
are  shown. 

The  AEH  approach  has  also  been  applied  to  biomechanics  applications.  Hollister  et  al. 
[56]  employ  the  homogenization  procedure  to  study  the  linear  elastic  material  properties  of 
porous  trabecular  bone  structures.  Global  and  local  information  are  estimated  and  apparent 
stiffnesses  were  computed  and  compared  with  experiments. 

As  an  extension  of  an  earlier  paper  for  the  damage  analysis  of  UD  fibers  [53],  Lene  and 
Paumelle  [57]  study  the  damaged  properties  of  woven  fabric  composites.  Again,  the  ap¬ 
proach  employs  the  elastic  constitutive  equations  with  the  damage  modeled  as  slipping  at 
an  interface  without  separation. 

The  homogenized  properties  from  the  AEH  approach  are  compared  with  those  from  me¬ 
chanics  of  materials  approaches  by  Hollister  and  Kikuchi  [58].  The  strain  energy  densities 
are  computed  and  compared  for  the  cases  when  the  homogenization  approach  is  employed 
for  a  coarse  FE  model  or  the  direct  solution  of  the  whole  refined  composite  is  studied.  The 
mechanics  of  material  approaches  are  based  on  flexural  and  stiffness  measurement  techniques 
that  represent  methods  of  estimating  the  bounding  material  properties  for  elastic  materials 
as  described  in  Chung  and  Tamma  [59]  which  provides  comparison  of  effective  properties 
for  other  types  of  homogenization  approaches  in  addition  to  AEH  and  standard  mechanics 
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methods.  Their  study  shows  the  agreement  between  AEH  and  other  approaches  for  linear 
elasticity. 


In  other  related  works,  Kawamoto  and  Kyoya  [60]  studied  porous  rock  and  pattern  bolting  of 
rocks  using  an  earlier  AEH  development  for  simple  elastic  materials  [55].  This  applications- 
oriented  investigation  computed  the  effective  elastic  material  properties  for  various  volume 
fractions  in  a  study  for  geomechanics.  Lefik  and  Schrefler  [61]  applied  AEH  to  linear  elastic 
beam  theory  including  shear  rotations  and  warping  using  FEM.  Equations  of  elasticity  for 
bricks  and  mortar  used  in  masonry  structures  were  homogenized  and  studied  by  Papa  [62]. 
The  effective  elastic  properties  and  the  damaged  properties  are  computed  via  homogeniza¬ 
tion.  Using  digital  image-based  (DIB)  geometric  modeling,  Golanski  et  al.  [63]  studied  the 
thermoelastic  behavior  of  layered  metal  matrix  heterogeneous  materials:  a  substrate  of  low 
alloy  steel  attached  to  a  composite  layer.  In  particular,  the  mismatch  in  material  properties 
between  successive  layers  and  the  microstresses  in  the  composite  layers  was  examined.  A 
review  of  the  method  for  linear  elasticity  is  presented  by  Chung  et  al.  [64]. 

The  elastic  problem  has  been  studied  extensively  in  the  literature.  Despite  the  seemingly 
simple  considerations  for  elastic  materials,  the  multidisciplinary  nature  of  sophisticated  prob¬ 
lems  and  the  complex  computational  issues  associated  with  the  homogenization  approach 
make  the  study  of  elastic  materials  a  cogent  beginning  for  more  sophisticated  engineering 
problems.  Some  of  these  are  described  next. 


2.5  Inelastic  Problems 


Stemming  from  early  mathematical  derivations  by  Sanchez-Palencia  [39]  and  Sanchez-Hubert 
[65],  numerous  efforts  have  investigated  the  curious  viscoelastic  behavior  of  heterogeneous 
materials  using  AEH.  Early  mathematical  literature  for  the  viscoelastic  problem  in  the  con- 
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text  of  AEH  shows  an  hereditary  effect,  the  dependence  of  the  present  state  of  deformation 
on  all  previous  states,  occurs  due  to  a  physical  and  mathematical  coupling  between  the 
non-hereditary  phases  and  the  elastic  and  viscous  material  tensors.  Phases  that  are  history 
independent  appear  to  produce  history  dependent  effects  when  composited  together.  Anal¬ 
ogous  behavior  in  other  rate-independent  inelastic  (such  as  elasto-plastic)  investigations  is 
not  observable  because,  regardless  of  the  length  scale  under  consideration,  the  incremental 
constitutive  (Prandtl-Reuss)  equations  are  Hookean  in  form.  That  is,  general  viscoelastic 
constitutive  laws  are  non-Hookean: 

P<x  =  Qe,  (43) 

where  P,  Q  are  differential  operators  whose  coefficients  are  material  properties  and  cr,  e 
are  stress  and  strain  tensors.  Other  types  of  inelastic  constitutive  models  can  typically  be 
expressed  in  a  Hookean  form: 

=  R(<r,  U)e,  (44) 

where  R  is  a  tensor  whose  terms  are  functions  of  material  properties,  stress  components 
(or),  internal  energy  ( U ),  or  other  state  variables.  A  Hooke’s  law  form  for  the  constitutive 
model  possesses  only  one  material  property  tensor  and  implies  that  the  AEH  approach  can 
be  extended  forthright  to  the  inelastic  problem.  The  homogenized  form  of  an  Hookean 
constitutive  model  is  Hookean  as  well.  Such  an  observation  is  not  permitted  in  non-Hookean 
forms.  As  a  result  of  the  equation  development,  mathematical  interactions  between  the 
multiple  terms  over  multiple  phases  are  likely  to  occur.  Another  interpretation  is  that 
physical  interactions  between  phases  occur  due  to  the  simultaneous  presence  of  solid  and 
fluid  effects. 

Sanchez-Hubert  and  Sanchez-Palencia  [65]  showed  mathematically  that  the  instantaneous 
viscoelastic  constitutive  equation  for  a  Kelvin- Voigt  material  can  be  related  to  long-term 
aging  or  long-memory  effects  via  the  AEH  formulation.  That  is,  the  Kelvin- Voigt  model 
applied  to  the  microlevel  constituents,  when  homogenized,  may  result  in  a  macrolevel  stress- 
strain  equation  different  in  form  from  that  of  the  original  microlevel  equation. 
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The  difference  is  a  long-term  aging  or  dissipative  corrector  term  that  manifests  itself  in  the 
macrolevel  equations.  Two  causes  mathematically  give  rise  to  the  corrector  term.  The  first  is 
the  asymptotic  expansion  of  the  displacement  variables,  the  central  assumption  of  AEH.  The 
second  is  the  presence  of  additive  rate-dependent  and  rate-independent  terms  in  the  stress- 
strain  relationship.  Sanchez-Palencia  [39]  shows  that  the  coupling  results  in  an  additional 
term  in  the  macro  stress  strain  relationship  whose  interpretation  is  that  of  long  memory. 
Unfortunately,  though  their  investigation  provides  insights  into  the  mathematical  origins  of 
the  long-memory  effects,  the  magnitude  of  those  effects  or  even  their  practical  demonstration 
were  omitted.  That  is,  the  relevant  equations  have  never  before  been  solved. 

Francfort  and  Suquet  [66]  later  showed  a  different  formulation  with  a  simpler  development  of 
the  viscoelastic  micro-macro  creep  equations.  Their  work  extended  the  earlier  developments 
[39]  by  investigating  convergence  to  the  continuum  solution  and  boundedness  of  the  principal 
quantities.  But  despite  the  depth  to  which  the  investigation  of  the  mathematics  is  presented 
in  Francfort  and  Suquet  [66]  and  the  novelty  of  the  formulations  in  Sanchez-Palencia  [39], 
many  of  the  computational  and  numerical  details  for  the  viscoelastic  multiscale  creep  problem 
in  practical  and  engineering  applications  have  only  recently  been  addressed  by  Chung  et 
al.  [67-69]. 

Viscoelastic  methods  via  AEH  for  heterogeneous  materials  have  been  attempted  and  are 
present  in  the  literature.  Two  efforts  are  available  [70,71],  described  below,  both  of  which 
make  several  assumptions  regarding  the  multiscale  issues. 

Dynamic  viscoelastic  response  was  investigated  where  the  effective  composite  properties  pre¬ 
dicted  were  compared  favorably  with  limited  measurements  [70].  The  dissipative  corrector, 
however,  is  avoided  by  imposing  an  appropriate  harmonic  oscillating  strain  comprised  of  time- 
dependent  and  time-independent  components  in  the  dynamic  viscoelastic  problem.  Through 
the  AEH  derivation,  only  the  time-independent  component  is  retained  hence  precluding  the 
time-dependent  dissipative  corrector  effect  from  appearing. 
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Shibuya  [71]  performed  the  asymptotic  expansion  of  the  displacements  in  the  Laplace  do¬ 
main  wherein  the  equilibrium  equations  governing  the  fiber  and  matrix  phases  are  assumed 
uncoupled.  The  governing  equations  are  rendered  in  the  Laplace  frequency  domain  where, 
together  with  AEH,  the  viscoelastic  response  in  a  single  time-step  approach  can  be  obtained 
for  Maxwell  materials  using  the  elastic-viscoelastic  correspondence  principle. 

Other  methods,  different  from  AEH,  that  treat  the  viscoelastic  problem  for  more  sophisti¬ 
cated  constitutive  models,  employ  analytical  methods  or  expeditious  assumptions  that  are 
inapplicable  to  general  composite  microstructures.  Thus,  despite  their  availability,  especially 
those  that  estimate  broad-spectrum  phase  behavior  to  which  the  quasi-static  problem  is  a 
subset,  the  important  limitation  is  in  their  assumption  for  simple  shapes  in  the  microstruc¬ 
ture.  An  additional  limitation  in  these  existing  methods  is  the  decisive  decoupling  of  the 
primary  length  scales.  The  AEH  method  still  remains  an  open  area  for  viscoelastic  materials. 

Research  in  the  homogenization  of  linear  problems  continues  to  be  of  interest  in  the  study 
of  heterogeneous  materials.  They  provide  the  mathematical  and  physical  foundations  upon 
which  extensions  to  more  complex,  nonlinear  problems  can  be  based.  Linear  problems  are 
also  important  in  computational  mechanics  because  important  implementational  and  com¬ 
puter  issues  are  easier  to  identify.  Ultimately,  however,  classes  of  realistic  multidisciplinary 
nonlinear  problems  are  of  greatest  interest  in  the  field  of  computational  mechanics.  Litera¬ 
ture  for  nonlinear  problems  employing  AEH  is  described  next. 


2.6  Nonlinear  Engineering  Problems 


Limited  AEH  efforts  are  in  the  literature  for  nonlinear  engineering  situations  due  to  the 
complicated  issues  associated  with  evolution  of  microlevel  properties.  To  date,  most  compu¬ 
tational  mechanics  efforts  to  treat  nonlinear  problems,  apart  from  the  present  study,  appear 
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exclusively  for  solid  mechanics.  Efforts  are  divided  between  geometric  and  material  nonlin¬ 
earity. 

In  perhaps  one  of  the  earliest  demonstrations  of  a  nonlinear  engineering  problem  using  com¬ 
putational  mechanics  and  AEH,  Guedes  [72]  employs  FEM  for  quasi-static  elasto-plasticity 
with  finite  deformations  but  without  updating  the  microlevel  properties  between  successive 
numerical  iterations.  The  true  elasto-plastic  problem  must  involve  change  of  material  prop¬ 
erties  in  the  microstructure  constituent  phases,  as  recognized  in  Suquet  [73].  As  yielding 
occurs,  the  material  properties  must  be  rehomogenized  to  accurately  model  the  evolving 
properties.  That  procedure  was  later  corrected  in  a  paper  by  Terada  and  Kikuchi  [31]. 

Ghosh  and  co-workers  [33,  34]  use  the  Voronoi  cell  Finite  Element  Method  (VCFEM)  to 
derive  the  governing  equations  for  two-scale  elasto-plastic  analysis  of  heterogeneous  materials 
with  AEH.  The  VCFEM  is  an  alternative  to  conventional  meshing  of  the  heterogeneous 
microstructure.  The  analysis  is  for  small  displacements  with  the  constitutive  model  expressed 
in  Hookean  form  with  an  instantaneous  tangent  modulus,  making  the  extension  from  AEH 
for  elasticity  straightforward.  The  computational  approach  employs  an  implicit  solution 
procedure  which  solves  the  nonlinear  governing  equations  iteratively. 

Galvanetto  et  al.  [74]  present  an  elasto-plastic  formulation  for  plane  problems  within  an 
asymptotic  expansion  framework.  By  attempting  to  treat  general  nonlinear  problems,  their 
investigation  poses  restrictions  on  the  forms  of  the  constitutive  equations  to  which  the 
method  can  be  extended.  Again,  it  is  apparent  that  constitutive  models  adhering  to  Hookean 
forms  can  account  for  heterogeneous  media  by  extending  the  AEH  approach  for  elastic  ma¬ 
terials.  This  approach  has  been  discussed  indirectly  in  other  investigations  [31]. 

Fish  et  al.  [32]  present  an  eigenstrain  formulation  to  account  for  Hookean  constitutive  laws 
with  additive  strains  in  an  elasto-plastic  framework.  A  two-point  integration  scheme  is 
proposed  to  decrease  computation  time  with  limited  loss  in  accuracy.  The  equations  are 
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integrated  implicitly  using  an  iterative  scheme. 


Other  work  in  nonlinear  AEH  are  varied.  Issues  of  geometric  nonlinearity  are  considered 
separately  by  Smit  et  al.  [75]  and  Takano  et  al.  [76].  Rate  dependent  nonlinear  composites 
have  been  considered  by  Wu  and  Ohno  [77].  And  dynamic  transient  problems  with  nonlinear 
elasto-plastic  models  have  been  reported  by  Chung  et  al.  [78]. 

Computational  methods  for  heterogeneous  media  via  AEH  are  still  in  a  state  of  early  devel¬ 
opment  as  evidenced  in  the  available  literature.  Many  areas  still  remain  as  fertile  ground 
upon  which  new  computational  methods  can  still  be  developed.  Several  conclusions  from 
this  literature  review  are  described  next. 


2.7  Demonstration  of  the  Homogenization  Method  in 
Elasticity 


Asymptotic  homogenization  approaches,  and  in  particular  their  associated  computational 
methods,  for  linear  problems  are  in  continued  development  as  evidenced  in  the  literature. 
Linear  problems  provide  the  sturdy  foundations  upon  which  nonlinear  phenomena  can  be 
studied. 

The  goal  of  this  section  is  to  present  AEH  in  a  simple  manner  for  the  elastic  problem. 
Elements  in  this  linear  formulation  provide  the  foundations  to  which  later  developments  will 
be  reference.  Also  discussed  are  the  FE  computational  issues  associated  with  the  model 
linear  elastic  problem.  These  will  ultimately  convey  to  any  extension  of  the  approach  for 
inelastic/nonlinear  problems  or  problems  outside  the  realm  of  continuum  scales. 

The  AEH  approach  for  linear  elasticity  is  well  established  and  rigorously  derived  in  the 
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literature.  This  section  provides  only  a  simple  derivation  to  demonstrate  the  approach  in 
practical  terms  for  computer  implementation.  In  section  2.7.1,  the  conventional  equations 
of  elasticity  are  defined  formally.  Then  in  section  2.7.2,  the  multiple  scale  equations  are 
derived.  In  section  2.7.3,  the  relevant  FE  formulations  are  described,  and  the  results  are 
presented  in  section  2.7.4. 

2.7.1  Equations  of  Conventional  Elasticity 

The  PDE  system  for  elasticity  is  described  here  for  a  bounded  connected  domain  Q,  €  5?3 
with  smooth  boundary  dQ  in  coordinates  (indices  denote  quantities  in  Jt3  and  the  Einstein 
summation  convention  is  used).  The  boundary  is  composed  of  two  parts,  and  <92Q,  which 
are  the  surfaces  of  non-zero  area  associated  with  the  Dirichlet  and  Neumann  conditions.  The 
body  £1  may  further  be  divisible  into  two  parts,  and  f 22,  by  a  smooth  surface,  T.  Figure 
9  shows  a  schematic  diagram  of  the  body  where  =  fli  U  fi2  and  dfl  =  dQ\  U  5f22. 


Figure  9.  Elastic  domain  definitions. 
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2. 7. 1.1  Boundary  Value  Problem  (BVP) 


The  BVP  is  comprised  of  four  key  components:  (1)  equilibrium,  (2)  compatibility,  (3)  consti¬ 
tution,  and  (4)  sufficient  boundary  conditions  for  well-posedness.  The  equilibrium  condition 
is  specified  by 

^+/S  =  0,  (45) 

where  is  the  Cauchy  stress  tensor  and  fj  is  the  body  force  vector.  In  the  framework  of 
linear  elasticity  with  displacements  as  the  unknowns,  compatibility  is  inherently  specified  by 
the  strain-displacement  relationship  given  by 


where  is  the  strain  tensor  and  Ui  is  the  displacement.  The  relationship  between  stress 
and  strain,  in  this  case  Hooke’s  Law,  gives  the  constitutive  equation  for  the  problem. 

Oij  —  Dijkitku  (47) 

where  DtJki  is  the  fourth  order  elasticity  tensor  which  is  piecewise  smooth  in  each  phase 
having  discontinuities  on  the  interface  T.  The  symmetry  and  ellipticity  of  are  specified 

by 

Dijki  —  Djiki  —  Dijik  —  Dknj ,  (48) 

Dijkitijtki  >  acijtij  ;oi>0  Ve^ (u)  (symmetric).  (49) 

Finally,  Dirichlet  and  Neumann  boundary  conditions,  corresponding  to  the  clamping  (uj) 
and  surface  force  traction  (Ft)  boundary  conditions  on  80.,  ensure  a  unique  solution.  These 
are  given  by 


U{  =  0 

(JijTlj  —  F) 
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on  9if2; 
on  82  £1. 


(50) 

(51) 


The  BVP  can  be  solved  in  one  of  many  number  of  different  approaches.  The  present  ap¬ 
proach  employs  FEM  to  approximate  the  solution  on  a  discretized  general  domain.  As  such, 
variational  forms  of  the  equations  expedites  the  formulations. 


2. 7. 1.2  Variational  Problem 

The  variational  formulation  of  the  elasticity  problem  is  presented  here  only  as  a  manner 
of  formality.  The  description  becomes  complete  when  presented  in  a  functional  analysis 
framework.  Additional  and  more  general  details  are  available  in  Bensoussan  et  al.  [29]  and 
Sanchez-Palencia  [39]. 

To  write  variational  forms  of  the  equations,  the  Hilbert  space  for  the  norm  of  (H1^))3  is 
defined  as 


V  =  {u;uie//1(fi)  ;ui|ftn  =  0}.  (52) 

The  symmetric  inner  product,  by  virtue  of  the  symmetry  of  the  material  property  tensor  in 
equation  (48),  is  given  by 

g(u,  v)  =  ^  Dijkieij(u)eki(v)dx  =  Dm^  (53) 

The  variational  form  of  the  BVP  given  in  equations  (45)-(47),  (50),  and  (51)  is:  Find  Ui  €  V 
such  that 

a(u,  v)  =  /  fiVidx  +  /  FiVids  Vvi  6  V.  (54) 

Jn  Ja2n 

The  variational  form  of  the  equations  can  then  be  extended  using  an  FE  approximation.  The 
FE  formulations  will  be  presented  in  the  next  section  in  the  context  of  general  heterogeneous 
media. 
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The  conventional  equations  of  elasticity  are  the  governing  equations  for  single  phase,  homo¬ 
geneous  materials.  Methods  for  heterogeneous  materials  typically  involve  an  homogenization 
step.  The  AEH  approach  involves  an  homogenization  step  that  is  derived  directly  from  the 
governing  equations.  This  characteristic  of  the  equations  and  their  developments  lends  math¬ 
ematical  rigor  to  the  homogenization  approach,  as  frequently  described  in  the  literature.  The 
“rigor”  will  be  described  in  detail  in  the  next  section. 


2.7.2  Homogenization  in  Elasticity 

This  section  describes  in  detail  the  approach  for  employing  a  perturbative  asymptotic  series 
in  e  for  the  displacement  variables  to  derive  multi-scale  equations  governing  elasticity.  We 
assume  that  the  material  under  consideration  has  a  microstructure  comprised  of  multiple 
phases,  and  the  orientation  and  shapes  of  the  phases  are  such  that  they  are  distributed 
repetitively,  periodically  in  all  three-dimensions  throughout  the  material. 

For  a  more  formal  mathematical  description,  let  Q  be  a  bounded  domain  in  5ft3  of  coordinates 
Xi  (or  x),  as  depicted  in  Figure  9.  The  domain  Q  denotes  the  macrolevel  or  global  structure. 
In  the  space  of  5R3  coordinates,  let  there  be  a  fixed  parallelepiped  Y  on  coordinates  yl  (or  y) 
of  edges  y°  (see  Figure  10).  All  other  parallelepipeds  in  the  body  are  obtained  by  an  integer 
translation  of  length  ny°  where  n  is  some  integer.  This  assemblage  of  parallelepipeds  implies 
that  a  single  parallelepiped,  or  unit  cell,  is  Y-periodic  in  y,  due  to  the  repetition  of  the  cell 
Y  throughout  the  body. 

In  a  conventional  homogeneous  material,  the  tensor  of  material  properties  is  indepen¬ 
dent  of  position  x.  If  the  material  is  not  homogeneous,  however,  the  material  property  will 
depend  on  x.  Moreover,  if  the  heterogeneous  material  contains  a  periodic  microstructure, 
where  the  length  of  the  period  is  much  smaller  than  any  other  lengths  appearing  in  the  prob¬ 
lem,  then  the  material  behavior  can  be  approximated  by  a  homogenized  material  property 
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Figure  10.  Y-periodic  parallelepiped  cell. 


tensor,  D^kl.  based  only  on  the  salient  features  of  one  period.  The  homogenization  method 
provides  a  means  of  computing  D^kl  such  that  the  global  problem  can  be  treated. 

The  periodicity  runs  throughout  the  entire  coordinate  system  The  microstructural  vari¬ 
ations  and  the  characterization  of  heterogeneous  information  are  defined  in  the  coordinate 
system  yt.  To  characterize  the  periodicity  of  the  material  properties,  e  is  defined  such  that 

£%«(*)  =  Dm(  5),  (55) 

where  e  is  a  real  positive  parameter  where  e  G  (0,  e0]-  The  function  D\^kl  is  eY-periodic  in 
Xi .  The  period  eY  is  a  repeated  array  of  parallelepipeds  in  3?3  whose  edges  are  length  £y°. 
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From  these  developments,  equations  (45),  (50),  (51),  and  (47)  are  rewritten  as 


dxe  ~~  ^  in 

(56) 

u\  —  0  on 

(57) 

oljTij  —  Fi  on 

(58) 

=  DtjkFkii  u£), 

(59) 

where  the  e  script  denotes  quantities  describing  the  true,  high-resolution  behavior  the  ma¬ 
terial  under  consideration.  The  given  forces  fi  and  iq  are  noticeably  not  periodic  because 
they  are  applied  on  the  surface  of  the  global  body  Q. 

To  summarize  thus  far,  quantities  such  as  the  stress,  cr?-,  with  the  script  e,  are  the  stresses 
in  x\  which  give  high-resolution  details  of  the  locally  periodic  nature  of  the  structure.  This 
stress  is  often  impractical  to  compute  due  to  the  large  amount  of  modeling  effort  needed  to 
create  a  full  size  mesh  with  all  the  details  of  the  microstructure.  Moreover,  the  computer 
time  needed  to  solve  the  problem  on  the  mesh  is  cost  prohibitive.  It  is  prudent,  therefore, 
to  introduce  an  homogenization  assumption.  This  is  done  by  identifying  the  global  body  in 
Xi  as  an  homogenized  body  with  a  definitive  unit  cell(s)  in  the  coordinate  system  y^.  The 
local  (or  micro)  coordinates  yi  are  related  to  the  global  (or  macro)  coordinates  aq  via  the 
parameter  e.  Thus,  the  coordinate  system  y,  is  a  magnified  scale  of  the  microstructure  whose 
absolute  lengths  are  given  in  the  x\  coordinate  system.  And  due  to  the  same  global  body 
defined  in  x\  and  X{,  the  distance  between  any  two  points  in  these  coordinate  systems  are 
equivalent. 

The  displacements  are  approximated  with  an  asymptotic  series  representation  in  e  given  by 

u\  (x)  =  (x,  y)  +  eu\l)  (x,  y)  +  e2  •  •  •  ,  (60) 

whose  terms  must  be  obtained  by  solving  sets  of  equations,  i.e.,  BVPs.  The  hierarchical 
equations  will  be  described  next.  The  Y-periodicity  of  the  material  requires  that  the  stresses, 
strains  and  displacements  also  be  Y-periodic.  The  functions  u^\x,  y)  (j=l,2,...)  are  assumed 
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smooth  functions  that  are  Y-periodic  in  the  variable  yi  such  that  yi  —  Xi/s  and  defined  for 
Xi  €  Q,yi  £  5ft3  independent  of  e.  Thus,  conversion  of  lengths  between  the  magnified  scale  y{ 
and  the  original  length  scale  x\  is  computed  using  e.  Henceforth,  it  is  common  practice  to 
refer  to  any  quantity  with  ^/-dependence  as  Y-periodic.  Any  quantity  that  is  solely  dependent 
on  x  with  the  script  e  notation  is  cY-periodic. 

Moreover,  derivatives  in  xf  must  now  be  expanded  in  a  chain  rule  for  the  space  x^  possessing 
a  periodic  microstructure  in  y{.  The  derivatives  are  now  over  two  length  scales  given  by 

A  =  A  +  Ii_.  (6i) 

dx\  dxi  e  dyi 


The  strains  in  equation  (46)  are  therefore  rewritten  as 


(u‘) 


dxi 


duf  1  du) 
+  - 


(0) 


+ 


duf'1 


(i) 


1  duf 


dxj  e  dyj  dxi  e  dy 
duf  duf  df 


dxj 


+ 


.(2) 


2  du\ 

+  £2-7rL-  +  e 


dVj 
duf 


+  S' 


dxi 


+ 


du 


(i) 


dyi 


dxj  dyj 


2  duf  duf 

+  £2^r^+£-  3 


dxi 


dyi 


+ 


_  i  ( .  M!' \.  i  (&?fL  +  ^  +  +  Ml) 

2e  ^  dyj  dyi  J  2  ^  dxj  dxt  dyj  dyt  J 

+  +  *£V... 

2  y  dxj  dxi  dyj  dyi  J 
=  je-71)(x,y)  +  €^)(x,y)+^1)(x,y)+e2---  . 


(62) 


46 


Hooke’s  Law  in  equation  (47)  also  can  be  written  in  the  expanded  form  as 


=  £«u(x)«iw  OO 


1  / Suf  ,  duf\^l(du?)duf  du?duf 


=  D,iU(y)  +  -g^j  +  2  (&T 

2  y  dxj  dxi  dxj  dtji 
=  -o\jl)  (x,  y)  +  cr|0)  (x,  y)  +  eg (x,  y)  +  e2  • 


J - £ - 1 - 

dyj  dyi 


Substituting  equation  (63)  into  (57)  gives  a  new  set  of  hierarchical  BVPs.  These  hierarchical 
BVPs  are  the  multiple  equations  relevant  to  each  scale  in  the  multiple  scale  problem  under 


consideration.  The  substitution  yields 


r-ia4~1}  -  E-2dA L. 

dxj  ^  dVj 


.odaii 


iggS!  |  £o?^L  |  £2^jL 

dxj  dy3  dx3 


,£i?±  +  ...  +  fi  =  o, 


and  upon  rearranging  terms,  yields 


.-2d<rH 


do^  +  d*l 

dxj  dyj 


do-- 

,  uuij 

dx3  dyj 


The  equations  must  be  valid  for  all  e  0.  Therefore,  the  coefficients  of  the  powers  of  e 
must  be  zero  identically.  This  implies  that  the  first  three  coefficients  are  written  as 


+  V^  =  0  ; 


dxj  dyj 


47 


(68) 


da. 


(o) 


da^ 


+  fi  ~  0 


dxj  dyj 

Equations  (66)-(68)  are  the  hierarchical  boundary  value  problems. 


In  summary  thus  far,  three  key  operations  have  been  performed.  In  the  first  operation,  the 
displacements  in  t-space  were  replaced  by  the  series  approximation  in  equation  (60).  In  the 
second  operation,  the  chain  rule  in  equation  (61)  was  used  for  derivatives.  The  chain  rule  is 
a  consequence  of  the  decomposition  of  the  true  structure  in  xf  into  an  homogenized  body  in 
Xi  (in  which  lengths  are  identical  to  those  in  the  coordinate  system  xf)  with  a  representative 
unit  cell  in  y*.  The  unit  cell  is  characterized  in  a  coordinate  system  that  is  a  magnified  scale 
of  the  microstructure  that  originally  resides  in  xf.  In  the  third  operation,  the  displacements 
and  chain  rule  were  substituted  into  the  conventional  BVP  for  elasticity,  and  by  setting 
the  coefficients  for  the  powers  of  e  to  zero,  three  equations  (66)-(68)  were  developed.  The 
first  and  second  equations  are  discussed  in  the  next  section  in  the  derivation  of  a  microlevel 
equation.  The  subsequent  section  employs  the  micro  equation  to  derive  a  coupled  multiscale 
BVP  for  linear  elasticity. 


2.7.2. 1  Micro  Equation 

The  micro  equation  refers  to  the  equation  governing  the  elastic  phenomena  on  the  domain 
characterized  by  the  ?/,  coordinate  system  (the  xz  coordinate  system  characterizes  the  macro 
domain).  The  micro  equation  is  derived  from  equations  (66)  and  (67)  where  the  first  equation 
defines  the  macro  displacement  variable  which  is  then  used  to  solve  the  second  equation 
to  give  the  first  order  perturbative  variable.  The  present  formulation  is  for  a  dual-scale 
heterogeneous  material  where  it  is  implied  that  continuum  level  descriptions  are  applicable 
to  both  length  scales. 

Substituting  for  the  first  three  orders  for  stresses  and  invoking  the  symmetry  associated  with 
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the  material  property  tensor  gives  the  following  alternate  forms  of  equations  (66)-(68). 


d  duf 

o  Dijkl  r, 

dyj  oyi 


=  0; 


(69) 


d  ^  du j;0)  d  ( duf  duf\ 

+  e^DiiU  (&T  +  A 


=  0; 


(70) 


An  ( ?^L_ 

dx3  m  \  dx, 


dujA  d  ( du ^  duf\ 


+  fi  =  0. 


(71) 


By  virtue  of  the  Y-periodicitv  of  uf\x,y)  and  the  ellipticity  (and  positivity  characteristics) 
of  Dijkii  it  is  evident  that 


duf ^ 

dyj 


=  0. 


(72) 


This  important  result  reveals  that  is  a  constant  in  yt,  or  stated  differently,  uf  is  the 
global  solution,  the  displacement  of  the  macrolevel  body  and  independent  of  y\.  With  this 
in  mind,  equation  (70)  is  simplified  to 


d  r\  ( duf 

Wi  m  A  +  ~w  j 


=  o. 


(73) 


Equation  (73)  is  an  equation  that  relates  the  zeroth  and  first  order  displacements.  Based 
on  the  global  solution,  the  perturbative  displacements,  can  be  obtained  from  the 
equation 


A  (n  duf  dDjjki 

dyj  y  zjkl  dyi  )  dxt  dyj 


(74) 


In  practice,  however,  the  solution  of  equation  (74)  is  cumbersome  in  an  FE  sense  if  uf 
varies  over  re*.  Such  a  scenario  requires  the  repeated  solution  for  uf  for  every  a  task 
that  is  computationally  onerous.  Or  stated  differently,  if  one  employs  a  straight-forward 
approach,  the  solution  of  equation  (74)  must  be  repeatedly  obtained  for  each  global  element. 
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Fortunately,  a  simpler  and  more  efficient  approach  is  available.  To  help  introduce  the  simpler 
version  of  the  same  problem,  the  Y-periodic  Hilbert  space  Vy  is  defined  as 

VY  =  {u  ;  Ui  €  H{oc(§ R3)  Y  -  periodic}  ,  (75) 

where  the  script  loc  refers  to  the  local  Y-periodicity  depicted  in  Figure  11.  In  Figure  11, 
the  solutions  at  points  a\  and  a2  are  approximately  equal  due  to  their  local  proximity  in 
Xj,  though  not  the  same,  and  their  respective  locations  within  each  period,  ?/*,  which  are 
the  same.  The  solutions  at  ai  and  /?2  are  different  because  their  locations  in  yi  are  entirely 
different  despite  their  local  proximity  in  The  solutions  at  /?2  and  ftz  are  also  different 
despite  their  similar  locations  yt.  The  difference  is  due  to  their  large  distance  apart,  much 
larger  than  the  size  of  a  single  period,  eyf. 


It  is  proposed  that  the  solution  of  the  perturbative  displacement  takes  the  form 


(i) 

u\  =  X; 


klduk 


(0) 


dxi 


+  tL(1) 


(x)> 


(76) 


where  uj^(x)  is  a  constant  of  integration  independent  of  y,  and  xf  is  the  solution  to  the 
auxiliary  variational  problem  given  by:  Find  xf  €  Vy  such  that 


r,  dxr  dVi 

Dijkl  ^  ^ 

oyi  oy-j 


dy  = 


f  Vj^p^iy 

Jy  dyi 


VVi  e  Vy. 


(77) 


The  function  xfl  is  often  termed  the  elastic  corrector.  Its  namesake  is  described  later  in  the 
development  of  the  homogenized  global  equations.  Other  terms  in  the  literature  describe  x\l 
as  the  characteristic  function  or,  more  classically,  its  gradient  is  the  so-called  misfit  strain. 


It  is  immediately  evident  that,  via  a  superposition  principle,  the  solution  of  equation  (77) 
eventually  yields  a  solution  to  equation  (74).  Equation  (74)  is  often  called  the  micro-equation 
because  it  relates  the  global  solution  to  the  perturbative  solution  and  is  solved  over  the 
domain  of  the  unit  cell.  It  is  also  the  equation  governing  the  smaller  of  the  two  length  scales 
and  illustrates  the  coupling  between  micro/local  («-^)  and  macro/global  (u-0"*). 
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Figure  11.  Local  periodicity. 


From  equations  (62)  and  (63),  approximate  microlevel  strains  and  stresses  can  be  computed 
using  the  solution  for  u The  expression  for  microlevel  strains  is  given  by  substituting 
equation  (76)  into  (62)  and  taking  the  zeroth  order  terms  of  e  which  gives 


4°W)  =  § 


1 

2 


1  /  du) 


=  u 


(0) 


dxj 


+ 


du) 


(0) 


dxj 


+ 


dxi 

du P 


duP  duf> 
+  -z2—  + 


dxj  dxi 


Hm^jn  + 


dxi 

dXr 


+ 


dxTn  du 


(0) 

m 


+ 


dxT  duf ' 
dxi 


dlLm 


dyj  J  dxn 


dyj  dxn  dyi 

(symmetric  about  ij). 


(78) 
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In  a  similar  way,  the  microlevel  stresses  can  by  obtained  from  the  expression 


40)(x>y)  =  Ai«(y)c«;(x,y) 


(0), 

tj 


—  Dijki(y)  |  SimSjn  + 


dxT 


dUm 


dt/j  )  &r„ 


(79) 


Localization  is  performed  by  employing  equation  (78)  to  yield  the  microlevel  strain  details. 
Similarly,  equation  (79)  gives  the  microlevel  stresses.  Localization  refers  to  the  process  of 
obtaining  microlevel  information,  or  information  on  the  y*  coordinate  system  from  the  global 
solution  (Wj-°^). 


In  summary  thus  far,  the  corrector  Xil  was  introduced  as  a  means  of  expressing  the  pertur¬ 
bative  displacement  in  a  closed- form  relationship  with  the  global  displacement  uf\  It 
is  of  key  importance  to  distinguish  between  the  global  displacement  uf'*  and  the  microstress 
aff.  Although  both  share  the  same  script  notation  for  the  zeroth  order,  the  gradient  of  the 
global  displacement  in  xf.  by  virtue  of  the  chain  rule  of  equation  (61),  does  not  immediately 
give  the  global  stress.  This  inequality  can  be  expressed  by 


duf'*  du. 


(o) 


dx] 


dxj 


(80) 


But  mathematical  consistency  of  the  formulations  can  be  demonstrated  by  showing  that  the 
gradient  of  in  x*  gives  the  volume  average  of  the  gradients  in  xf.  This  is  characterized 
by 


where  |F|  is  the  volume  of  the 
(79)  gives  the  global  stress. 


-L^<>(x,yMy  =  |b,  (81) 

unit  cell.  That  is,  the  average  of  the  microstress  in  equation 
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2. 7. 2. 2  Single-Point  Localization  Paradox 

Equations  (78)  and  (79)  exemplify  a  key  deficiency  of  the  AEH  approach.  Despite  the 
equations  being  “micro  equations,”  length  scale  information,  namely  the  scale  parameter 
£,  is  conspicuously  absent.  This  precludes  the  exchange  of  length  scale  information  across 
multiple  scales.  The  absence  of  e  is  attributable  to  the  first  order  approximation  of  the 
asymptotic  series  in  e.  The  gradients  of  the  first  order  displacement  for  strains  and  stresses 
are  therefore  zeroth  order  in  s.  The  physical  interpretation  is  that  the  AEH  approach  involves 
point-wise  localization.  Upon  homogenization,  the  scale  and  details  of  the  microstructure 
are  “smeared”  away.  When  one  returns  to  compute  microlevel  information  via  localization, 
the  microstrains  and  stresses  are  computed  at  a  single  point  in  the  global  body.  An  entire 
unit  cell  exists  at  the  single  point,  x *.  Thus,  the  point-valued  global  strain  is  then  used  to 
compute  the  local  strains  (and  stresses)  over  all  t/j. 

A  paradox  is  therefore  evident.  The  scale  parameter  e  is  a  quantity  that  magnifies  the 
coordinate  system  of  the  microstructure.  But  the  microstructure  is  assumed  mathematically 
disconnected  from  the  global  body.  The  scale  parameter  is  not  used  to  relate  a  position 
in  yi  to  the  original  heterogeneous  body  in  the  e-space.  Thus,  if  two  points  (x|  and  xf) 
are  selected  arbitrarily  close,  both  will  have  an  entire  unit  cell  Y.  The  two  unit  cells  may 
therefore  overlap,  hence  the  paradox.  We  presently  call  this  the  single-point  localization 
paradox  and  is  depicted  schematically  in  Figure  12. 

It  is  reassuring,  however,  that  in  equation  (78),  the  strain  term  outside  of  the  parentheses 
is  the  only  term  dependent  on  x*.  The  other  terms  are  constants  or  depend  on  y*.  Thus,  in 
the  limit  as  the  two  points  approach  each  other  (x)  — >  xf),  the  global  strain  at  the  point  is 
single- valued.  Even  if  two  selected  points  in  x*  are  very  close  and  have  two  different  unit  cells 
associated  with  them,  the  microstrains  computed  over  those  two  unit  cells  become  the  same 
as  x\  -»  xf.  That  is,  if  x|  and  xf  are  very  close  but  not  the  same,  the  difference  in  computed 
microstrains  may  be  trivially  small.  An  analogous  argument  can  be  made  for  microstresses. 
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Figure  12.  The  single-point  localization  paradox.  The  approximations  of  the 
AEH  approach  remove  length  scale  information  at  the  microlevel. 
This  allows  for  unit  cells  to  overlap  conceptually  and  leads  to  the 
localization  paradox. 


The  single-point  localization  paradox  implies  that  equations  (78)  and  (79)  are  only  quan¬ 
titative  “estimates”  of  the  local  strain  and  stress.  To  have  an  exact  prediction  entails  the 
solution  of  the  fully  refined  problem  with  all  microscale  details  of  the  structure  modeled  us¬ 
ing  a  prohibitively  large  number  of  degrees  of  freedom.  Such  is  the  trade-off  between  solving 
an  homogenized  solution  with  estimates  of  the  local  behavior  and  an  extremely  large  single 
scale  high  resolution  problem  possessing  a  large  number  of  degrees  of  freedom.  The  local 
estimates  are  merely  representative  of  the  behavior.  In  an  FE  context,  this  local  region  is 
an  element  or  the  region  near  the  quadrature  point  where  the  gradient  terms  of  equations 
(78)  and  (79)  are  computed. 


2. 7.2. 3  Macro  Equation 

Equation  (76)  is  now  used  to  complete  the  derivation  of  the  effectively  homogenized  equations 
for  linear  elasticity.  That  is,  the  effective  macro  BVP  is  formulated  here.  Substituting 
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equation  (76)  into  (71)  gives 


+  q —  Dijkl 

dy3 


du^j  duf} 

dxi  dyi 


(82) 


+  /»  —  o. 


If  up"1  is  to  be  Y-periodic,  then  by  virtue  of  the  definition  of  Y-periodicity,  equation  (82) 
admits  a  unique  solution  for  vjp  up  to  an  additive  constant  if  and  only  if 

jy  ■kDi>t‘a-^dy  =  “•  (83) 

Substituting  equation  (82)  into  (83)  yields 


du^  dx™n  dum 
dxt  dyi  dxn 


dy  =  0 


(84) 


and  upon  rearrangement  and  dividing  both  sides  by  the  volume  of  the  unit  cell  |Y|  gives 

^{w\Ld'iu 

Equation  (85)  and  equations  (50)  and  (51)  constitute  the  homogenized  BVP  in  elasticity. 
The  material  property  tensor  is  now 


dyn, 

SkmSin+  Xk 


dyi 


dy>-t+f<=° 


(85) 


Dijmn  ~  jy  |  ^  Dijkl  ^km^ln  + 


gxT 
dyi  J 


dy. 


(86) 


Thus,  the  average  zeroth  order  stress,  also  the  macrolevel  effective  stress,  is  given  by 


—  Hh 

—  Uijkl 

—  Dh 
~  U ijkl 


duf^ 

~d^' 


where  the  volume  average  notation  is  implied  by  the  brackets: 


(87) 

(88) 


Equation  (86)  provides  the  motivation  for  referring  to  xf  as  a  corrector  function  because 
of  its  corrective  feature  in  the  equation.  In  its  absence,  equation  (86)  gives  the  standard 
volume  average  of  the  material  properties,  characteristic  of  the  Rule-of-Mixtures  approach 
for  estimating  homogenized  properties  of  heterogeneous  materials.  Thus  the  corrector  effec¬ 
tively  “corrects”  for  the  presence  of  multiple  phases  within  the  unit  cell  and  their  collective 
interactions  everywhere  else  in  the  global  body  via  periodic  boundary  conditions  on  dY . 
The  gradient  of  xf  can  also  be  interpreted  as  a  misfit  strain  because  it  accounts  for  the  mis¬ 
match  strains  resulting  in  the  multiple  phase  microstructure.  However,  this  study  chooses 
the  “corrector”  terminology  for  its  physical  interpretation  where  xf  corrects  the  heteroge¬ 
neous  material  from  the  otherwise  conventional  homogeneous  understanding  of  the  governing 
equations. 

This  concludes  the  derivation  of  the  multiscale  BVP  for  linear  elasticity.  The  equation 
for  the  microscale  is  given  in  equation  (77)  and  the  equation  for  the  macroscale  is  given 
in  equation  (85).  Solving  equation  (77)  gives  xf  which  is  then  used  in  equation  (86)  to 
calculate  the  effective  homogenized  properties.  The  homogenized  properties  are  used  in 
equation  (85)  to  compute  the  global  solution,  Once  the  global  solution  is  obtained, 

the  perturbative  displacement  term,  uf  \  can  be  computed  using  equation  (76).  Moreover, 
the  microlevel  information  can  be  computed  using  xf  and  the  localization  equations  given 
in  equation  (78)  for  strains  and  equation  (79)  for  stresses.  The  homogenization-localization 
cycle  demonstrated  provides  the  integral  link  allowing  causal  relation  of  the  continuum  BVP 
across  multiple-length  scales. 

Although  the  mathematical  formulations  are  well-defined  and  discussed  extensively  in  the 
literature,  the  practical  implementational  issues  are  often  nebulous.  In  the  next  section,  FE 
equations  are  formed  and  presented  along  with  a  computational,  implementational  procedure 
which  provides  the  conceptual  understanding  required  to  employ  AEH  in  general  numerical 
methods. 
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2.7.3  Finite  Elements,  Elasticity,  and  the  Homogenization  Ap¬ 
proach 

This  section  provides  a  practical  approach  for  computing  the  corrector,  the  effective 
material  properties,  D^kl,  and  the  localized  perturbative  displacement  term,  iq-  ,  in  the 
framework  of  FEM.  The  description  and  the  utility  of  the  AEH  approach  differ  subtly  from 
the  earlier  mathematical  descriptions  when  placed  in  the  context  of  numerical  methods.  This 
section  attempts  to  discuss  these  differences  and  present  a  simple  approach  for  employing 
the  method  for  relatively  complicated  heterogeneous  problems. 


2.7.3. 1  Computing  the  Corrector 


The  solution  of  the  variational  equation  given  in  equation  (77)  gives  the  corrector,  xf •  The 
variational  equation  is  rewritten  here  for  convenience:  Find  Xil  €  Vy  such  that 


L 


D, 


ijkl ' 


dx™71  dvi 
dyi  dy-j 


dy 


=Lvi 


dDijmn 


\fvt  e  Vy. 


For  simplicity,  consider  the  material  under  investigation  to  be  in  the  class  of  orthotropic 
materials  with  nine  independent  components  in  the  material  property  tensor.  Henceforth, 
Voigt  notation  will  be  employed  which  exploits  the  symmetry  of  the  strains  and  allows  the 
following  notational  simplification  to  be  made: 
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Furthermore,  the  appropriate  piecewise  linear  Sobolev  spaces  are  implied  for  the  test  and 
trial  functions.  A  less  rigorous  engineering  description  is  favored  here  over  a  mathematical 
FE  procedure  in  terms  of  function  spaces.  Additional  details  can  be  found  in  numerous  FE 
references  [55,79-82].  The  standard  FE  representation  for  element  strain,  {e},  is  given  by 

W  =  [£]{«},  (91) 

where  [£]  is  the  element  strain  matrix  and  {u}  is  the  vector  of  nodal  displacements.  The 
element  stresses  are  given  by 


{a}  =  [D][B}{u}, 


(92) 


where  [D]  is  the  matrix  of  material  properties  and  {a}  is  the  vector  form  of  the  symmetric 
stress  tensor.  Thus,  the  FE  analogue  of  equation  (77)  is  given  in  element  form  by 

f  [B]T[D]{B]dye[x]  =  f  [B]T[D)dy’,  (93) 

JYe  JY* 

where  the  script  e  denotes  element  quantities  associated  with  the  discretized  FE  domain  of 
the  unit  cell,  namely  the  body  Ye.  Equation  (93)  can  also  be  written  for  a  3-D  problem  as 


m  w  =  [fd] 

["Sot  x  ndof)  ("dot  x  61  Kof  x  6) 

where 

[Ke]  =  [  [BflDPJdy*  and  [FD]  =  [  [B\T[D]iy‘, 

JYe  JYe 

and  n[jof  is  the  number  of  nodal  degrees  of  freedom  of  each  element  in  Ye. 


(94) 


(95) 


2. 7.3. 2  The  RHS  Matrix  [FD] 

The  corrector  is  noticeably  a  matrix,  not  a  vector  as  is  the  case  for  displacements  in  conven¬ 
tional  elasticity.  This  is  a  result  of  the  multiple  right-hand  side  vectors  signifying  the  load, 
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[, Fd ].  Each  column  of  the  load  “matrix”  corresponds  to  a  column  in  the  material  property 
matrix,  [£>].  Stated  differently,  equation  (93)  is  actually  a  set  of  six  matrix  equations  for  3-D 
elasticity.  These  can  be  summarized  in  a  conventional  sense  of  solving  for  a  single  column 
of  unknowns  x ^  as 


\Kt]{xin)}=  [  [B\T{D^}dy\ 

Jy * 

where  n  =  1  •  •  •  6,  implying  a  total  of  six  solutions  of  the  problem,  with 
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The  six  solutions  provide  the  six  columns  of  [x]  in  equation  (94). 


Thus,  for  a  given  elastic  problem  of  a  heterogeneous  body,  the  approach  requires  the  solution 
of  six  sets  of  equations  to  compute  the  corrector.  Each  solution  corresponds  to  a  particular 
“mode”  or  “character”  of  the  unit  cell:  three  normal  modes  and  three  shear  modes. 

To  obtain  a  unique  numerical  solution,  a  zero  constraint  on  x  must  be  applied  somewhere 
within  the  FE  model  of  the  unit  cell.  This  can  be  done  arbitrarily  because  only  the  corrector 
gradient  is  needed  in  obtaining  the  homogenized  properties.  For  explicit  computation  of  the 
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perturbative  displacement  ti>  ,  it  is  optimal  to  choose  a  node  exploiting  the  symmetry  of  the 
parallelepiped.  In  the  present  study,  the  zero  constraint  is  applied  on  a  corner  node  where 
any  three  edges  intersect. 


2. 7.3. 3  Periodic  Boundary  Conditions 

Additionally,  boundary  conditions  must  be  applied  to  simulate  the  periodicity  as  required 
by  the  AEH  approach.  Periodicity  is  enforced  by  requiring  nodal  value  equality  on  opposing 
boundaries  of  the  unit  cell  model.  To  illustrate  this,  consider  the  3-D  parallelepiped  as 
depicted  in  Figure  13.  Boundary  conditions  for  all  degrees  of  freedom  at  the  relevant  nodes 
on  opposing  faces  of  the  following  form  must  be  applied 

x*(o»  ife,  2/3)  =  Xi(y°,  V2, 2/3),  (98) 

x*(yi>  o,  yz)  =  Xi(yu  2/2^3),  (99) 

Xi(yi,y2,o)  =  Xi(yi>y2,y|)-  (10°) 

In  practice,  a  simple  linear  algebra  operation  on  the  assembled  stiffness  matrix,  [AT],  is  all 
that  is  required  to  impose  these  boundary  conditions.  Hinton  et  al.  [83]  provide  an  approach 
that  reduces  the  order  of  the  stiffness  matrix  in  accomplishing  this  task. 

The  periodicity  assumption  is  crucial  to  the  numerical  procedure  because  it  imposes  the  con¬ 
dition  originally  specified  theoretically  in  equations  (55)  and  (75).  The  explicit  dependence 
of  the  equations  on  and/or  y*  are  predicated  on  the  definition  of  X  in  3?3  and  the  periodic¬ 
ity  definition  of  Y.  These  conditions  are  tantamount  to  enforcing  the  assumption  that  each 
given  unit  cell  is  embedded  in  an  array  of  identical  unit  cells.  Methods  of  enforcing  other 
types  of  periodicity  conditions  are  discussed  in  Whitcomb  [26]  and  Chung  [84].  Though  no 
studies  in  the  literature  appear  to  investigate  the  merits  of  one  set  of  periodicity  conditions 
over  another,  it  traditionally  understood  that  the  periodicity  breaks  down  near  boundaries 
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Xi(  a  )  =  Xi(a’) 

X,(b)  =  Zi(b’)  i-U.3 
X,(c  )  =  X,(c’) 

Figure  13.  Periodic  boundary  conditions.  Nodal  values  for  x  on  opposing 
faces  are  equal. 

or  regions  of  global  discontinuities  [37]. 

However,  the  assumption  of  periodicity  is  necessary  in  sophisticated  approaches  because  it 
indirectly  accounts  for  the  presence  and  interactions  of  the  other  unit  cells.  Without  such 
conditions,  only  dilute  concentrations  of  inclusions  are  permissible,  the  correctness  of  the 
present  formulations  notwithstanding. 

Finally,  to  obtain  a  solution  with  rigid  body  modes  removed,  a  zero  displacement  condition 
must  be  applied  arbitrarily  in  the  body  Ye.  The  effective  properties  are  independent  of  the 
location  where  the  zero  condition  is  applied.  The  relative  magnitudes  of  the  nodal  values  of 
[x],  however,  are  easiest  to  interpret  when  the  zero  condition  is  applied  at  a  location  that 
exploits  the  cubic  symmetry  of  the  unit  cell  in  Figure  13.  The  present  study  enforces  this 
condition  at  an  arbitrary  corner  where  any  three  edges  intersect.  Then,  by  virtue  of  the 
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periodicity  assumption,  all  corners  also  have  zero  displacement. 


Once  the  corrector,  [x],  is  determined  at  the  nodes  of  the  unit  cell  FE  model,  the  effec¬ 
tive  property  matrix  of  the  model  can  be  computed.  In  the  next  section,  an  approach  for 
computing  the  homogenized  effective  elastic  properties  of  the  unit  cell  is  described. 


2.7.3.4  Computing  the  Homogenized  Property  Matrix 

The  computation  of  the  homogenized  property  matrix  in  FEM  comes  from  a  direct  modifi¬ 
cation  of  equation  (86)  which  gives 

^elm  t r 

[BA]  =  E?r-[  o']([/]  +  [B']M),  (wi) 

e=l  Ktot 

where  [. Dh ]  is  the  homogenized  property  matrix  for  the  unit  cell,  neim  is  the  number  of 
elements  in  the  unit  cell,  Ve  is  the  volume  of  the  element,  Vtot  is  the  total  volume  of  the 
unit  cell,  [/]  is  the  identity  matrix,  [Be]  is  the  element  strain  matrix,  and  [xe]  is  the  matrix 
of  nodal  [x]  values  connected  to  that  element.  Summing  over  all  the  elements  provides  the 
volume  average  over  the  unit  cell  with  the  added  correction  that  [xe]  provides.  The  script  e 
in  all  of  the  terms  denotes  the  micro  element  in  Ye. 

The  last  term  of  equation  (101)  is  the  so-called  corrector  term  which  accounts  for  the  presence 
of  the  inhomogeneities  in  the  microstructure.  It  also  accounts  for  the  interactions  between 
a  given  cell,  Ye,  and  the  periodic  array  of  cells  surrounding  it.  This  is  what  is  referred 
to  earlier  in  section  2.7.2. 1  as  the  misfit  strain,  in  keeping  with  micromechanical  elasticity 
theories  [85,86]. 

The  gradients  of  [xe]  are  computed  at  the  integration  points  of  the  elements  in  Ye  where 
they  are  numerically  most  accurate.  For  linear  problems,  the  gradients  at  the  integration 
points  can  be  averaged  to  provide  a  single  gradient  for  the  entire  element.  The  gradient  of 
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[xe]  in  equation  (101)  is  this  average  in  the  present  context  of  elasticity. 

It  is  of  interest  to  note  that  the  integral  in  equation  (86),  due  to  the  division  by  the  unit 
cell  volume,  adds  no  dimensionality  to  the  homogenized  properties.  This  implies  that  the 
effective  property  matrix  is  only  a  function  of  the  relative  lengths  of  inclusions  and  matrix. 

The  effective  material  property  matrix  is  then  employed  in  the  global  equations  as  the 
material  property  matrix  of  the  global  element.  It  is  therefore  implied  that  each  global 
element  possesses  one  unit  cell  representative  of  the  microstructure  in  that  global  region  of  the 
structure.  After  the  global  solution  is  obtained,  the  local  gradients  within  the  representative 
unit  cell  can  be  computed.  The  next  section  describes  the  approach  for  computing  this 
“localization”  information. 


2. 7. 3. 5  Computing  the  Perturbative  Term 

The  perturbative  term,  is  determined  at  the  nodes  of  the  body  in  Ye.  The  term  is 
a  function  of  the  global  deformation  gradient  which  varies  from  one  global  element  to  an¬ 
other.  Thus,  if  the  gradients  through  the  global  body  are  not  uniform,  then  the  perturbative 
displacements  will  also  vary  accordingly.  This  can  also  be  interpreted  as  the  local  periodic¬ 
ity  assumption  permitting  variations  of  the  microlevel  information  from  one  location  in  the 
global  body  to  another. 

In  the  integrated  multiple  scale  approaches  of  this  study,  the  local  information  is  computed 
using  the  deformation  gradient  at  the  numerical  integration  points  of  the  global  element. 
Thus,  the  integration  points  of  the  global  elements  serve  as  conduits  of  information  through 
which  passes  the  macro  information  to  the  micro  domain. 
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In  equation  (76),  the  perturbative  displacement  was  given  as 

,(i) 


u) 


X 


ki^u 

dx 


(o) 

— +41,(x)- 


The  boundary  conditions  used  to  compute  xf  described  in  section  2. 7.3.3  implies  that  the 
constant  of  integration  of  equation  (76),  is  zero.  Then,  the  nodal  perturbative  displace¬ 
ments  in  a  3-D  problem  are  given  by 


{«<■>}  =  ix)  m 

{wdof  x  1}  [^dof  x  6]  [6  x  n^of] 

where  {it^}  is  the  vector  of  nodal  perturbative  displacements  at  the  microlevel,  [x]  is  the 
matrix  of  nodal  corrector  values  at  the  microlevel,  [B9]  is  the  strain  matrix  of  the  global 
element  for  which  the  microlevel  unit  cell  is  being  considered,  and  tidof  total  number 

of  degrees  of  freedom  in  the  FE  mesh  for  the  body  in  Ye. 


{<,) 


(102) 


Finally,  the  microlevel  strains  and  stresses  can  be  computed  using  equations  (78)  and  (79) 
from  the  perturbative  displacements  given  in  equation  (102).  In  the  next  section,  the  specific 
procedural  steps  are  outlined. 


2. 7.3.6  Computational  Procedure 

The  procedural  steps  for  homogenization  and  localization  of  a  linear  elastic  problem  are 
described  here.  It  is  assumed  for  simplicity  that  the  global  geometry  possesses  only  one 
microstructural  representative  unit  cell.  That  is,  the  FE  problem  involves  the  solutions  of 
the  governing  equations  over  two  meshes  -  the  global  body  mesh  and  the  local  unit  cell  mesh. 
Should  the  microstructure  vary  at  different  locations  within  the  global  body,  it  is  understood 
that  the  procedure  must  be  extended  to  account  for  multiple  unit  cell  meshes. 

In  the  following  steps,  the  dimensions  of  the  matrices  are  provided  for  a  3-D  scenario.  The 
matrix  sizes  can  be  changed  according  to  the  dimensional  of  the  problem.  The  distinction 
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between  global  FE  quantities  and  local  quantities  are  specified  by  script  g  or  e,  respectively. 


Step  1.  Create  the  FE  model  of  the  global  body  in  Xe. 

Step  2.  Create  the  FE  model  of  the  microstructural  body  in  Ye. 

Step  3.  Solve  for  the  corrector  in  the  variational  equation  given  by  equation  (77)  or  the 
discretized  corrector  in  equation  (93): 

jy.  mT  m  mdy*  [X]=jye  mT  mdY* 

Kof  x  6]  [6  X  6]  [6  X  ned0{]  [nedoi  x  6]  [nedoi  x  6]  [6  x  6] 

The  RHS  terms  are  given  in  equation  (97)  and  the  symmetric  boundary  conditions 
are  imposed  according  to  equations  (98)-(100). 


Step  4.  Compute  the  effective  elastic  tensor  by  substituting  the  correctors  determined  in 
step  3  into  equation  (86)  or  (101): 

P‘]  =  E2r5fc  P1(  W+  [B'l  M) 

[6  x  6]  [6  x  6]  [6  x  6]  [6  x  n‘dJ  fa*  x  6]  ' 

The  result  of  this  step  is  the  homogenized  material  property  tensor  or  homogeniza¬ 
tion. 

Step  5.  Solve  the  global  problem  for  using  the  homogenized  properties  of  step  4  subject 
to  boundary  conditions  and  loads  on  the  global  body.  The  solution  is  obtained  from 
equation  (85)  or  the  discretized  form  given  by 

/*«  mT  m  m**  wr  (103) 

(ndof  x  6]  [6  x  6]  [6  X  ned0{]  {n^of}  Kof  x  3]  {3} 

The  values  of  from  this  step  give  the  traditional  global  solution  analogous  to  the 
conventional  response  of  an  homogeneous  material. 
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Step  6.  Compute  the  average  local  strains  at  the  desired  discrete  locations  in  Xe  (usually  at 
quadrature  points)  via  equation  (78)  on  an  element-by-element  basis  over  Ye  which 
is  given  by 

{e<°> }  (x,  y)  =  ([/]  +  m  M)  [B’\  {«<»> }  ,  (104) 

where  the  script  g  and  e  are  used  to  denote  the  FE  strain  matrices  for  the  global 
element  and  the  local  element,  respectively.  Note  that  the  gradient  term,  du^/dx 
or  [B9]  {«(«)},  is  constant  in  y?-.  This  means  that  for  a  given  microstructural  unit 
cell,  [j B 9)  is  a  constant. 

Step  7.  Compute  the  local  stresses  via  equation  (79)  also  on  an  element-by-element  basis 
over  Ye : 

{<r<°>}  (x,y)  =  [!>']  ([/]  +  [Be]  M)  [B»]  {«<»>}  .  (105) 

The  material  property  tensor  [D]  in  equation  (105)  varies  according  to  y*.  This  step, 
together  with  step  6,  completes  the  localization  procedure  that  yields  the  local  strain 
and  stress  information. 

2.7.4  Comparisons  With  Other  Homogenization  Methods 

This  section  is  divided  into  two  subsections.  The  first  subsection  discusses  the  similarities 
between  AEH  and  classical  homogenization  approaches  to  develop  an  analogy  to  facilitate 
the  understanding  of  the  AEH  homogenization  approach.  In  the  second  subsection,  sim¬ 
ple  comparisons  between  the  present  homogenization  approach  and  other  homogenization 
approaches  are  presented  in  the  context  of  inplane  plate  stiffness. 


2. 7.4.1  Classical  Approach  Analogy 


There  exists  a  rich  body  of  literature  in  the  understanding  of  effective  properties  of  com¬ 
posites  [13,85-89].  Many  so-called  classical  approaches  employ  analytical  techniques  which 
place  assumptions  on  the  simply  shaped  geometry  of  the  microstructure.  The  assumption  on 
geometry  and  isotropy  of  the  phases  permits  the  development  of  closed  form  approximations 
for  the  effective  properties.  They  consequently  preclude  considerations  for  complex  substrate 
geometries,  such  as  woven  fabrics,  and  are  limited  to  spherical,  cylindrical,  ellipsoidal  or  other 
types  of  simply  shaped  inclusions.  Others  depart  from  a  classical  methodology,  and  instead, 
opt  for  numerical  techniques  or  analytical  methods  derived  from  classical  approaches  to  ac¬ 
commodate  the  complex  microstructure  of  advanced  high  performance  composites  [24-28,38]. 
Some  approaches  ostensibly  claim  causal  links  between  the  local  and  global  behaviors  [38] 
much  like  the  advocated  AEH  approach  of  this  study. 

However,  there  are  three  key  distinctions  that  set  the  present  AEH  approach  apart  from 
other  computational  homogenization  methods.  The  first  is  the  ability  of  the  AEH  approach 
to  provide  homogenized  properties  for  linear  and  nonlinear  continuum  BVPs.  The  second  is 
the  capability  of  estimating  local  response  at  very  small  length  scales  within  the  framework 
of  seamlessly  integrated  multiscale  mathematical  derivations.  Approaches  such  as  Woo  and 
Whitcomb  [38]  provide  a  means  of  estimating  global  and  local  response  but  use  kinematic 
balances  between  length  scales  based  on  no  mathematical  development  that  inherently  re¬ 
stricts  the  schemes  to  the  elastic  regime.  The  third  is  the  similarities  between  the  AEH 
approach  and  classical  homogenization  approaches.  Analogies  with  classical  approaches  can 
be  drawn  due  to  the  rigorous  equation  development  of  the  AEH  approach,  and  this  is  the 
topic  of  study  in  this  section.  Although  it  is  beyond  the  scope  of  the  present  study  to  review 
all  classical  homogenization  approaches  -  elasticity,  heat  transfer,  etc.  -  there  are  inherent 
features  over  a  broad  range  of  homogenization  schemes  that  are  of  interest  for  comparison 
purposes  with  the  present  computational  method. 
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The  comparisons  and  analogies  drawn  here  are  only  to  qualitatively  examine  the  correctors 
described  in  earlier  sections  in  a  classical  mechanics  context.  The  equation  morphology  of 
correctors  can  draw  parallels  within  the  framework  of  the  general  theory  of  eigenstrains  [85] 
and,  for  a  more  specific  case,  Eshelby’s  formulation  [86,90].  The  rigorous  treatment  of  these 
parallel  observations  is  the  subject  of  another  investigation.  The  objective  here  is  to  loosely 
show  the  physical  inferences  permitted  by  these  comparisons  for  didactic  purposes. 


2. 7.4. 1.1  General  Theory  of  Eigenstrains  In  a  strict  sense,  an  eigenstrain  refers  to 
any  general  non-elastic  strain  such  as  thermal  expansion,  phase  transformation,  or  misfit 
strains  [85].  Presently,  in  the  context  of  the  correctors  [x]  which  is  employed  to  effectively 
“correct”  the  homogenized  material  properties  by 

^elm  t  r 

it  is  evident  that  computing  the  average  stresses  based  on  the  effective  elastic  properties 
yields 

^elm  ry 

{*(0)}  =  £  (M  +  [seM  -te(0)>-  (107) 

e=l  ^tot 

In  Mura’s  general  theory  of  eigenstrains,  the  fundamental  assumption  regarding  linear  elas¬ 
ticity  is  that  the  total  strain  (e)  is  the  linear  superposition  of  the  elastic  strain  (ee)  and  the 
eigenstrain  (e*),  given  by  (for  notational  simplicity,  the  zeroth  order  script  of  the  asymptotic 
series  is  henceforth  omitted) 

{e}  =  {ee}  +  {e*}.  (108) 

The  elastic  strain  is  related  to  the  stress  by  Hooke’s  Law  which  gives 

{a}  =  [£)]{£«}  =  [D]  ({e}  -  {€-}) .  (109) 

By  comparing  equation  (107)  with  equation  (109),  it  is  immediately  evident  that  the  second 
strain  component  upon  distributed  multiplication  in  equation  (107)  is  analogous  with  the 
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eigenstress  component  in  equation  (109).  The  eigenstress  is  the  constitutive  counterpart  of 
the  aforementioned  eigenstrain.  In  the  present  context  of  heterogeneous  materials,  the  eigen- 
strain  is  precisely  the  misfit  strain  which  accounts  for  the  inhomogeneous  stress  distribution 
that  arises  from  the  presence  of  an  inclusion  in  a  homogeneous  matrix.  As  a  consequence, 
in  keeping  with  the  physical  descriptions  of  Mura  [85],  the  stress  and  strain  fields  are  “dis¬ 
turbed”  from  the  trivial  state  in  the  AEH  approach  due  to  the  fictitious  eigenstrain  (re: 
corrector)  in  the  homogeneous  material  that  accounts  for  the  additional  inclusion  phase  and 
the  subsequent  interactions  between  other  inclusions. 


2.7.4. 1.2  Eshelby’s  Formulation  A  more  specific  case  of  the  eigenstrain  interpretation 
is  seen  by  parallel  comparison  with  Eshelby’s  formulation  for  isotropic  inclusions  [86,90].  In 
Eshelby’s  original  formulations,  the  inclusions  and  matrix  are  of  the  same  material.  The 
premise  of  the  approach  is  to  equate  the  strain  energy  of  a  body  in  which  a  local  region  is 
subject  to  inhomogeneous  loads  (inclusion)  with  that  of  an  homogeneous  body  containing  a 
particular  distribution  of  body  forces  to  account  for  the  inclusion.  Misfit  strains  are  those 
that  result  from  the  body  forces  that  account  for  the  presence  of  the  inclusion,  thereby 
disturbing  the  homogeneous  matrix  from  its  initial  state  and  inducing  a  strain  field. 

The  problem  of  an  inhomogeneous  material  loaded  externally  by  an  applied  field  is  decom¬ 
posed  into  the  sum  of  two  problems.  The  first  is  that  of  a  homogeneous  material  subject  to 
the  same  external  loads  as  the  original  problem.  The  second  is  the  homogeneous  material  in 
the  absence  of  external  fields  subject  to  a  misfit  strain  (e*)  that  accounts  for  the  presence  of 
the  inclusion.  The  induced  strain  (e)  caused  by  the  presence  of  the  inclusion  can  be  obtained 
by  transforming  the  misfit  strain  field  via  the  transformation  given  by 

{«}  =  [£]{<*}.  (no) 

where  [E]  is  the  so-called  Eshelby  tensor.  The  decomposition  of  the  problem  implies  that 
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the  internal  energies  are  additive.  By  comparing  the  energies  of  the  various  problems,  the 
true  energy  of  the  heterogeneous  problem  was  derived  as  the  sum  of  the  internal  energy  of 
the  same  problem  with  the  inclusion  removed  plus  an  additional  interaction  energy.  That  is, 

U  =  U0  +  if/iNT,  (HI) 

where  U  is  the  energy  of  the  heterogeneous  problem,  UQ  is  for  the  homogeneous  problem 
where  the  matrix  has  been  removed,  and  UmT  is  the  interaction  energy.  The  interaction 
energy  for  a  problem  with  applied  tractions  is  defined  by 

UlnT=  [(tS-ttuDdS,  (112) 

where  S  is  the  surface  of  the  inclusion  and  U  and  u{  are  the  surface  tractions  and  the 
deformation,  respectively.  The  script  0  denotes  the  fictitious  problem  where  only  the  matrix 
is  considered  in  the  absence  of  the  inclusion. 

In  the  AEH  approach,  in  the  context  of  FEM,  the  total  internal  energy  is  given  as  the  sum 
over  Ye  of  the  internal  energy  in  each  element  given  by 

u=J2  '« {<<0)}T p]  (M  +  MVl)  {f(0)} .  (113) 

at  elm 

where  the  ellipticity  and  symmetry  of  the  terms  have  been  used.  In  a  more  abbreviated 
form,  equation  (113)  can  be  written  as 

U  =  UD  +  Ucorr,  (114) 

where  the  energy  components  are  self-evident  by  comparison  to  equation  (113).  Via  non¬ 
rigorous  manipulations,  it  can  be  shown  that  equation  (114)  becomes 

U  =  (UD  +  Um)  -  (Um  -  UCorr)  ,  (115) 

where  Um  is  a  “catch-all”  term  whose  value  can  be  tuned  by  comparing  with  equation  (111) 
such  that  U0  =  UD  +  Um.  Then,  it  is  clear  that  the  interpretation  of  the  misfit  strain  is 
applicable  to  the  AEH  approach. 
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Clearly  the  derivation  for  the  corrector  of  the  present  AEH  approach  is  different  from  that 
for  the  interaction  energy  in  Eshelby’s  formulations.  But  the  additional  interaction  energy 
shows  precise  parallels  in  the  interpretation  of  the  equation  forms  and  lends  credence  to  the 
“corrector”  interpretation  of  x  in  equation  (77).  This  superficial  comparison  of  the  equations 
is  meant  only  to  reveal  the  additive  nature  of  the  corrector  and  show  its  consistency  with 
interpretations  of  classical  homogenization  approaches. 

Moreover,  other  computational  approaches  such  as  those  defined,  for  instance,  by  Foye  [24, 
25],  Whitcomb  [26],  Gowayed  and  Yi  [27],  and  Woo  and  Whitcomb  [38],  are  not  as  easily 
comparable  to  classical  approaches.  Though  the  mechanics  of  material  schemes  employed 
therein  are  more  intuitive  from  a  force-balance  perspective,  it  is  precisely  for  this  reason  they 
cannot  be  extended  to  the  nonlinear  regime  and  show  the  rigorous  coupling  across  length 
scales. 

In  the  next  section,  quantitative  comparisons  are  undertaken  to  show  the  agreement  of  the 
AEH  approach  with  other  numerical  and  analytical  homogenization  approaches  for  plate¬ 
like  composites  in  linear  elasticity.  Though  the  AEH  approach  possess  features  not  found  in 
other  computational  homogenization  approaches,  it  is  still  expected  that  for  linear  problems 
agreements  with  existing  results  can  be  obtained. 


2. 7.4. 2  Numerical  Comparisons 

In  this  section,  quantitative  comparisons  are  presented  between  several  homogenization 
methods.  The  homogenized  inplane  stiffness  of  a  woven  fabric  plate  structure  is  used  for 
the  comparison.  The  homogenized  results  from  the  AEH  approach  are  compared  with  the 
so-called  Strain  Energy  Balance,  Finite  Element  Plate  Approximation,  and  Area  Averaging 
approaches.  In  addition,  two  analytical  methods  due  to  Ishikawa  and  Chou  [91]  for  the 
so-called  mosaic  and  crimp  models  are  used  to  further  verify  the  results  with  well-accepted 
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theories  for  woven  fabric  composites.  The  purpose  here  is  to  demonstrate  the  different 
effective  elastic  properties  of  composite  materials  obtainable  using  other  homogenization 
methods  and  to  show  the  definitive  localization  feature  of  the  AEH  approach.  The  distinc¬ 
tion  between  the  present  work  and  other  homogenization  approaches,  therefore,  is  illustrated 
quantitatively.  The  results  presented  in  this  section  are  directly  from  a  published  paper  [59] 
and  the  associated  thesis  [84]  on  the  effective  homogenized  elastic  properties  of  woven  fabric 
composites. 

The  present  considerations  are  for  general  composites  and  therefore  focus  on  methods  which 
can  homogenize,  say,  woven/undulating  stiff  cylindrical  inclusions  embedded  in  a  soft  matrix. 
In  light  of  the  scope  of  this  section,  homogenization  approaches  which  can  estimate  effective 
properties  of  complex  geometries  are  considered.  In  principle,  these  approaches  are  similar 
to  classical  mechanics  approaches  such  as  those  described  in  Refs.  [85, 86]  in  that  only  ho¬ 
mogenization  is  feasible.  This  illustrates  the  key  component  of  AEH  absent  in  other  existing 
approaches  -  AEH  possesses  the  ability  to  establish  a  causal  link  between  the  continuum 
BVPs  across  the  relevant  length  scales. 

A  full  comparative  study  of  selected  homogenization  schemes  was  conducted  using  the  finite 
element  meshes  of  the  microstructure  depicted  in  Figures  14(a)-(j).  The  study  compares 
only  the  inplane  Young’s  modulus  because  of  the  limitations  inherent  in  using  the  plate 
approximation  (see  Chung  [84]).  Figures  15(a)— (b)  show  the  moduli  for  the  mosaic  and 
crimp  models  and  compare  the  results  from  the  four  homogenization  schemes.  The  mosaic 
and  crimp  models  are  theoretical  idealizations  of  the  woven  geometry  of  fibers  used  in  woven 
fabric  composites.  They  were  originally  proposed  by  Chou  [92]  for  analytical  methods. 

Values  for  ng  used  were  2,  3,  4,  5,  10,  and  oo.  It  is  clear  by  observation  of  these  geometries 
that  they  are  implicitly  50%  fill  tows  by  volume  and  50%  warp  tows.  The  longitudinal  (fill) 
and  transverse  (warp)  material  properties  of  the  fiber  bundles  are  representative  of  an  epoxy 
graphite  composite  in  Table  1. 
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(a)  Mosaic  model,  ng  =  2 


(b)  Crimp  model,  ng= 2 


(c)  Mosaic  model,  ng= 3 


(d)  Crimp  model,  ng= 3 


(e)  Mosaic  model,  ng= 4 


({)  Crimp  model,  ng=  4 


(g)  Mosaic  model,  ng—b 


(h)  Crimp  model,  ng= 5 


(i)  Mosaic  model,  ng= 10  (j)  Crimp  model,  np=10 

Figure  14.  Finite  element  meshes  of  mosaic  and  crimp  models  for  various 
weaviness  parameter,  ng. 

The  results  indicate  that  the  AEH  approach  for  linear  elasticity  yields  results  within  the 
bounds  of  the  analytical  mosaic  and  crimp  models.  Furthermore,  the  AEH  numerical  re¬ 
sults  agree  well  with  the  other  computational  homogenization  approaches.  These  results  are 
indicative  of  the  good  agreement  expected  for  other  classes  of  linear  problems.  Nonlinear 
problems  typically  can  be  decomposed  into  a  series  of  smaller  linear  problems  in  numerical 
implementation.  In  light  of  this,  the  mathematical  consistency  of  the  formulations  further 
suggest  the  viability  of  AEH  for  treating  problems  in  the  nonlinear  regime. 


(a)  Mosaic  Model 

-  Crimp  Model  Theory 


(b)  Crimp  Model 

Figure  15.  Comparisons  of  inplane  moduli  for  the  mosaic  and  crimp  models 
for  varying  weaviness  parameter,  1  jng. 
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Table  1.  Constituent  elastic  material  properties  for  epoxy/graphite 
composite. 


Property 

Longitudinal 

Transverse 

Ex  (GPa) 

113.0 

8.82 

E,  (GPa) 

8.82 

8.82 

Ez  (GPa) 

8.82 

8.82 

vxy 

0.0235 

0.495 

Vyz 

0.495 

0.495 

VXz 

0.0235 

0.495 

Gxy  (GPa) 

4.46 

2.95 

Gyz  (GPa) 

2.95 

2.95 

Gxz  (GPa) 

4.46 

2.95 

Consider  now  the  global  problem  depicted  schematically  in  Figure  16.  The  FE  model  of 
this  geometry  is  depicted  in  Figure  17  with  2,541  nodes  and  2,000  elements.  Two  types  of 
microstructures  are  employed  in  this  section.  The  first  microstructure  is  representative  of 
an  orthogonal  non- woven  fiber  composite  depicted  in  Figure  18  with  1,917  nodes  and  1,472 
element.  The  second  model  is  representative  of  a  plain  woven  fabric  composite  depicted  in 
Figure  19  with  1,819  nodes  and  1,408  elements.  Constituent  material  properties  for  this 
example  are  listed  in  Table  2.  Each  material  component  is  assumed  isotropic.  The  objective 
here  is  to  compute  both  the  global  deformation  and  stresses  while  simultaneously  computing 
the  approximate  local  stresses  at  the  element  locations  depicted  in  Figure  16.  The  stresses 
are  approximate  due  to  the  homogenization  assumption  which  smears  the  inhomogeneous 
medium,  simulating  homogeneity,  and  the  subsequent  localization  which  uses  the  solutions 
obtained  from  the  smeared  behavior.  Element  no.  1,110  is  located  in  a  region  that  is 
expected  to  be  of  lower  overall  stresses  than  element  no.  1,990.  The  bar  will  experience 
stress  gradients  due  to  the  three-dimensionality  of  the  problem.  The  dimensions  of  the  bar 
are  2mxlmxlm. 
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Figure  18.  FE  model  of  (micro)  orthogonal  non-woven  fiber  composite. 


(a)  2D  plain  weave  composite 


(b)  Fiber  component  (c)  Matrix  component 

Figure  19.  Plain  weave  composite  (micro)  FE  models. 
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Table  2.  Constituent  properties. 


E  (GPa) 

G  (GPa) 

V 

Epoxy  resin 

3.5 

1.3 

0.35 

E-Glass 

72.0 

27.7 

0.30 

Ti-metal/?21-S 

112.0 

41.8 

0.34 

SCS-6 

393.0 

157.2 

0.25 

For  the  3-D  orthogonal  non-woven  and  the  2-D  plain  weave  composites,  the  elastic  material 
tensors  are  computed.  For  the  orthogonal  non-woven  fiber  composite,  the  elastic  tensor  is 
given  by 

21.1  5.3  5.3  0  0  0 

5.3  21.1  5.3  0  0  0 

5.3  5.3  21.1  0  0  0 

0  0  0  3.4  0  0 

0  0  0  0  3.4  0 

0  0  0  0  0  3.4 

For  the  2-D  plain  weave  composite,  the  elastic  tensor  is  given  by 

28.5  5.0  8.3  0  0  0 

5.0  12.3  5.0  0  0  0 

8.3  5.0  28.5  0  0  0 

0  0  0  3.3  0  0 

0  0  0  0  3.3  0 

0  0  0  0  0  6.9 

The  corrector  provides  a  measure  of  the  heterogeneity  of  the  microstructure,  namely,  the 
interaction  between  the  inclusion  and  the  matrix.  In  Figures  20  and  21  the  “mode”  shapes 
are  illustrated  via  isocontours  where  the  grey-scales  represent  the  magnitudes  of  the  corrector 
term. 
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(a)  Mode  11 


(b)  Fiber  cut-out  in  mode  11 


(c)  Mode  22 


(d)  Fiber  cut-out  in  mode  22 


(e)  Mode  33 


(f)  Fiber  cut-out  in  mode  33 


Figure  20.  Corrector  mode  shapes  for  orthogonal  non- woven  microstructure 
(magnified  10%  of  model  scale). 
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(c)  Mode  23 


(d)  Fiber  cut-out  in  mode  23 


(e)  Mode  13 


(f)  Fiber  cut-out  in  mode  13 


Figure  20.  Corrector  mode  shapes  for  orthogonal  non-woven  microstructure 
(magnified  109c  of  model  scale) (continued). 
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(e)  Mode  33  (f)  Fiber  cut-out  in  mode  33 

Figure  21.  Corrector  mode  shapes  for  plain  weave  microstructure  (magnified 
10%  of  model  scale). 
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(e)  Mode  13 


(f)  Fiber  cut-out  in  mode  13 


Figure  21.  Corrector  mode  shapes  for  plain  weave  microstructure  (magnified 
10%  of  model  scale) (continued). 


In  Figure  20(b)  the  fibers  extend  in  the  ^-direction  corresponding  to  the  x11  mode  shape. 
The  fibers  expand  in  the  11-direction  as  a  result  of  its  stiffer  behavior.  This  is  in  accordance 
with  equation  (77)  where  the  “loads”  are  defined  by  the  magnitude  of  the  property  difference 
across  the  phase  boundaries.  In  this  case,  the  stiff  fiber  has  higher  elastic  constants  which 
“loads”  the  softer  matrix  in  the  11-direction.  In  the  event  where  the  matrix  is  of  a  stiffer 
elastic  material  than  the  inclusion,  the  difference  in  properties  is  negative  indicating  that 
the  matrix  “loads”  the  fibers.  In  either  event,  the  loading  gives  the  appearance  of  the  stiffer 
material  expanding  while  the  softer  material  contracts.  Analogous  observations  can  be  made 
for  Figures  20(d)  and  (f),  and  Figures  21(b),  (d),  and  (f). 

Noteworthy  in  Figures  20  and  21  are  the  periodic  shapes  of  the  modes.  Each  unit  cell  can 
be  placed  adjacent  to  identical  cells  in  all  three  dimensions.  This  characteristic  reflects  the 
Y-periodicity,  or  the  translational  symmetry  of  the  unit  cell.  They  appear  to  fit  seamlessly 
together,  as  in  a  jigsaw  puzzle,  making  a  continuous  array.  The  symmetry  was  enforced 
mathematically  by  application  of  the  periodic  boundary  conditions. 

The  global  bar  twisting  solutions  are  illustrated  in  Figures  22(a)  and  (b)  for  the  non- woven 
and  woven  microstructures,  respectively.  A  line  load  of  P  =  10  kN/m  is  applied  counter¬ 
clockwise  around  the  edges  of  the  bar’s  tip.  The  associated  Von  Mises  stresses  and  strains 
are  illustrated  in  Figures  23(a)  and  (b)  and  Figures  24(a)  and  (b).  The  isotropic  effective 
properties  of  the  non-woven  microstructure  yields  isotropic  macro  behavior  of  the  bar.  The 
plain  weave  microstructure  is  transversely  isotropic,  and  this  behavior  is  reflected  in  Figures 
23(b)  and  24(b)  at  the  tip  where  the  respective  stresses  and  strains  are  not  distributed 
uniformly  around  the  perimeter  of  the  cross  section.  The  stresses  and  strains  at  the  bar  tip 
with  the  non-woven  microstructure,  however,  are  distributed  evenly. 

In  addition  to  the  standard  global  behavior,  which  is  obtainable  using  conventional  elasticity 
theory,  the  present  approach  provides  a  means  of  estimating  the  microstress  response.  The 
microstresses  are  computed  at  the  two  locations  indicated  in  Figure  16.  The  first  location  at 
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Figure  22.  Final  bar  deformation  shapes  and  magnitude  contours  (in  mil¬ 
limeters).  g4 
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Figure  23.  Von  Mises  stress  magnitudes  (in  Pascals). 
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(b)  Plain  weave  microstructure 

Figure  24.  Yon  Mises  strain  magnitudes  (xl0~3m/m). 
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element  no.  1,110  is  near  the  center  axis  of  the  bar  which  is  the  region  of  lowest  stress  in  the 
cross-section  plane.  The  second  location  at  element  no.  1,990  is  closer  to  the  outer  surface 
of  the  bar  where  it  is  expected  to  experience  larger  stresses.  Although  the  microstresses 
can  be  approximated  anywhere,  the  two  locations  are  chosen  presently  to  demonstrate  the 
differences  in  microstresses  using  the  AEH  approach.  In  contrast,  a  similar  problem  involving 
a  high  resolution  mesh  would  require  significant  modeling  effort  and  solution  time. 

The  microstress  contours  for  the  appropriate  geometries  are  depicted  in  Figures  25  and  26. 
The  differences  in  the  stresses  are  evident  between  Figures  25(a)  and  (c)  in  the  legends 
to  the  right  of  each  figure.  The  stresses  in  Figure  25(c)  are  clearly  greater  than  those  in 
Figure  25(a).  Whereas  the  median  microstress  at  element  no.  1,110  is  8.115  Pa,  the  median 
microstress  at  element  no.  1.990  is  53.9  Pa.  Analogous  observations  can  be  made  for  their 
respective  fiber  cut-outs  in  Figures  25(b)  and  (d)  or  for  the  plain  weave  microstructure  in 
Figure  26. 

The  average  of  the  microstresses  over  the  unit  cell  gives  the  global  stress  because  the  Y- 
periodic  correctors  cancel  when  volume  averaged  in  equation  (79).  This  is  what  gives  the 
AEH  approach  mathematical  consistency  between  scales. 

Finally,  although  the  present  discussions  have  not  fully  exploited  the  microstructural  infor¬ 
mation  obtained  by  the  present  approach,  further  studies  can  be  conducted  to  interpolate 
greater  refinement  of  microlevel  detail.  Figure  27  depicts  a  possible  area  of  investigation  to 
extract  stress  details  at  smaller  scales  by  identifying  and  manipulating  stress  information 
inside  the  cell.  This  topic  can  be  treated  in  the  context  of  scientific  visualization  in  future 
research. 

In  contrast  to  a  the  linear  elastic  problems  considered  in  this  section  which  involves  a  sin¬ 
gle  step  of  homogenization  and/or  localization,  time-dependent  problems  may  require  re¬ 
peated  homogenization.  This  is  especially  in  cases  where  the  material  properties  of  the 
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Figure  25.  Approximate  microstresses  for  non-woven  microstructure  (in  Pas¬ 
cals)  . 

constituents  change  with  time.  In  the  next  section,  a  linear  form  of  a  time-dependent  prob¬ 
lem  is  considered  in  the  viscoelastic  regime  where  the  BVP  is  a  function  of  time  due  to  the 
time-dependence  of  the  constitutive  equation.  However,  despite  the  unchanging  nature  of 
the  linear  material  properties,  the  global  problem  still  requires  the  rehomogenization  of  a 
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Figure  26.  Approximate  microstresses  for  plain  weave  microstructure  (in 
Pascals). 

limited  number  of  quantities  due  to  non-intuitive  effects  that  occur  when  time-dependent 
multiple  term  constitutive  models  are  homogenized. 
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Figure  27.  Additional  microstress  details  by  further  refinement  and  inspec¬ 
tion  of  the  microstructure. 

2.8  Conclusions  From  the  Literature  Review 


Four  primary  conclusions  can  be  drawn  from  the  literature  for  AEH:  (1)  computational 
methods  for  linear  and  elastic  engineering  problems  are  well  established  and  based  on  math¬ 
ematically  consistent  derivations  found  in  the  mathematics  literature;  (2)  some  problems 
originally  posed  and  discussed  cursorily  in  the  mathematics  literature,  such  as  in  Bensous- 
san  et  al.  [29]  and  Sanchez-Palencia  [39],  have  yet  to  be  fully  investigated,  understood,  and 
demonstrated/proven  for  practical  situations;  (3)  computational  methods  for  inelastic  and 
nonlinear  mechanics  remain  fertile  areas  for  growth  and  novel  developments,  and  (4)  despite 
the  ostensible  label  of  multiscale.  AEH  is  traditionally  applied  only  to  scales  where  con¬ 
tinuum  equations  apply,  and  no  efforts  show  applicability  to  scales  smaller  such  as  at  the 
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atomistic  level. 


It  is  clear  that  AEH  remains  an  active  area  of  research  for  computational  mechanics.  In  the 
next  section,  the  AEH  method  serves  as  the  framework  for  a  new  method  of  linking  atomistic 
scales  to  continuum  scales.  To  date,  no  documented  efforts  have  shown  either  theoretically 
or  in  application  that  this  is  possible. 
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3.  Formulation  of  a  Multiscale 
Atomistic-Continuum  Homogenization 
Method 


3.1  Overview 


A  large  amount  of  interest  has  recently  focused  on  the  multiscale  problem  involving  atoms 
and  continua.  It  is  widely  accepted  that  many  effects  on  the  continuum  germinate  at  the 
atomic  level.  Events  such  as  fracture,  fatigue,  and  inelastic  material  response  can  be  traced 
back  to  the  evolution  of  the  atomic  structure  of  the  material. 

Methodologies  for  linking  a  continuum  to  an  atomistic  domain  can  be  found  in  the  literature 
as  early  1971  [1].  Finite  element  methods  were  later  employed  in  Mullins  and  Dokainish  [2] 
using  a  numerically  decoupled  domain  approach  with  spatially  overlapping  atomistic  and 
continuum  regions.  A  review  of  some  of  these  methods  can  be  found  by  Cleri  et  al.  [3].  Among 
these  early  analytic  and  computational  studies,  frequent  issues  regarding  the  treatment  of  the 
interface  arose  which  were  primarily  handled  through  creative  use  of  kinematic  constraints. 

More  recently  Tadmor  et  al.  [4]  developed  an  entirely  FE-based  formulation,  the  so-called 
quasicontinuum  method.  Similar  efforts  were  made  through  the  so-called  handshaking  or 
coupling-of-length-scales  (CLS)  method  by  Broughton  et  al.  [5]  by  increasing  the  atomic 
resolution  via  the  tight-binding  (TB)  method.  And  the  dynamic  problem  was  studied  with 
a  generalized  scaling  approach  in  coarse-grained  molecular  dynamics  (CGMD)  by  Rudd  and 
Broughton  [6]  to  better  handle  the  propagation  of  waves  through  the  atomistic-FE  interface 
and  the  FE  far  field. 
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Multiscale  methods,  such  as  those  previously  described,  have  traditionally  been  limited 
mainly  to  localized  regions  of  interest.  For  example,  the  applications  to  which  these  methods 
have  been  applied  involve  small  sets  of  dislocations  and  cracks  and  very  few  limited  analy¬ 
ses  of  their  mutual  interactions.  The  localized  regions  on  which  these  simulations  are  run 
typically  span,  at  most,  several  microns.  The  limiting  assumption  in  these  works  is  the  use 
of  kinematic  constraints  to  tie  together  the  equations  and  disparate  length  scales.  Driving 
the  resolution  of  the  discretized  continuum  finite  elements  intrinsically  restricts  the  size  of 
the  continuum  and  leads  to  smaller  overall  dimensions  of  the  problem  which  can  only  be 
overcome  by  large  use  of  computer  resources.  Furthermore,  the  kinematic  constraints  on 
the  atoms  and  continua  lead  to  incompatibility  issues  arising  at  the  interface  such  as  ghost 
forces  [7]. 

The  asymptotic  expansion  homogenization  method  has  been  widely  studied  by  applied  math¬ 
ematicians  for  many  years.  Numerous  authoritative  texts  on  the  basic  theory  can  be  found  in 
the  literature,  for  instance,  by  Bensoussan  et  al.  [29],  Sanchez- Palencia  [39],  and  Bakhvalov 
and  Panasenko  [93].  And  despite  the  prolific  research  in  the  field,  no  attempts  have  been 
documented  for  extending  the  theory  to  atoms. 

In  this  part  of  the  report,  a  computational  framework  for  homogenization  of  the  atomistic 
problem  is  presented.  Using  two  concurrent  domains,  one  for  the  macroscale  continuum 
domain  and  one  for  the  atomic  scale  domain,  concurrent  self-consistent  sets  of  equations  are 
derived.  Atoms  in  arbitrary  configurations  and  structures  of  unlimited  size  are  permitted. 
Through  the  asymptotic  expansion  homogenization  technique,  a  set  of  hierarchical  equations 
are  derived  based  on  hyperelasticity.  At  the  local  level,  the  atomistic  equations  are  used  under 
the  assumption  of  the  harmonic  approximation  to  generate  the  effective  properties  needed 
to  solve  the  effective  global  level  equations.  The  Cauchy-Born  rule  [12]  is  applied  to  the 
atoms  to  enforce  the  gross  deformation  of  the  continuum  on  the  atoms.  This  circumvents 
the  need  to  apply  kinematic  constraints  by  making  use  of  the  weak  averaging  properties  of 
homogenization  and  causes  the  small  scale  equations  to  be  linear. 
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The  contents  of  this  part  are  as  follows.  In  section  3.2,  the  conventional  continuum  equations 
are  shown  eventually  leading  to  a  variational  form  based  on  the  principle  of  virtual  work. 
Then,  in  section  3.3,  the  multiscale  equations  are  developed  resulting  in  two  sets  of  equations 
which  govern  the  local  and  global  length  scales.  By  introducing  the  atomistic  potential  in 
sections  3.4  and  3.5  the  details  of  the  atomistic  formulations  are  presented  and  are  cast 
in  a  variational  form  for  use  in  the  multiscale  homogenization  method.  In  section  3.6,  the 
derivatives  of  the  atomistic  energy  potential  needed  to  complete  the  derivation  of  the  method 
are  provided  in  a  general  form.  In  section  3.7,  1-D  demonstrative  examples  are  worked  and 
shown.  Closing  remarks  are  discussed  in  section  3.8.  Additional  details  of  the  derivative  of 
the  Tersoff-Brenner  type  II  potential  are  shown  in  the  Appendix. 


3.2  Continuum  Formulations 


This  section  describes  the  kinematics,  stress  definitions,  and  linear  momentum  conservation 
laws  needed  to  develop  the  homogenization  method  from  atomistic  principles. 


3.2.1  Kinematics 


Consider  an  open  set  V  in  9?3  that  deforms  to  the  configuration  v  in  5ft3.  Points  in  V  are 
denoted  X  =  (A"i,  A'2,  A3)  €  V  and  are  called  material  points,  while  points  in  v  are  denoted 
x  =  (x\,X2,  X3)  €  v  and  are  called  spatial  points.  The  deformation  is  a  one-to-one  mapping 
through  (j>  so  that  x  =  <f>( X).  The  deformation  gradient  is  defined  by 


_d(pdx 
F  ~  dX~  dX~  °X’ 


F  ■  = 

13  dX. 


d(j)i  dxi 


dXj' 


(118) 


where  Vo  signifies  the  gradient  taken  with  respect  to  V.  The  determinant  of  F  is  termed 
the  Jacobian  and  is  defined  by  J  =  detF.  The  right  Cauchy-Green  strain  tensor  is  defined 
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by 


C  =  FtF,  (119) 

and  the  Green  strain  tensor  is  defined  by 

E  =  j(C-I).  (120) 


3.2.2  Stress  and  Equilibrium 


The  material  representation  for  the  conservation  of  linear  momentum  is  defined  by 

Vo  •  P  +  f0  =  0,  (121) 

where  P  is  the  first  Piola-Kirchoff  stress  tensor  and  f0  is  the  body  force  per  unit  of  undeformed 
volume.  In  rate  form,  it  is  given  by 

Vo  •  P  +  f0  =  0.  (122) 


Using  the  principle  of  virtual  work,  equation  (122)  can  be  rewritten  as 

/  (v0  •  pj  S-adV  +  J  f0  ■  6udV  =  0,  W5u, 


(123) 


where  <5u  is  the  virtual  displacement.  Then,  using  the  definition  for  traction  with  respect  to 
the  undeformed  body,  equation  (123)  can  be  rewritten  as 


/  P:Vo<5udF=  /  t0  •  SudA  +  /  f0  •  6udV. 

JV  JdV  JV 


(124) 


We  invoke  the  notion  of  hyperelasticity  by  assuming  that  the  atomistic  potential,  W,  which 
is  a  function  of  the  atom  positions,  can  be  expressed  in  terms  of  strain.  Intrinsically,  this 
assumes  that  the  strain  energy  density  (or  the  free  energy  at  zero  temperature)  is  equivalent 
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to  the  atomistic  energy  potential.  Following  classical  continuum  mechanics,  one  can  then 
define  the  first  Piola-Kirchoff  stress  as 


P  = 


dW 
dF  ’ 


P.  _ 
ij  ”  dFv' 


and  the  first  Lagrangian  elasticity  tensor  as 

c_  d2w  dp 


Ciikl  — 


d2W  dPij 


(125) 


(126) 


dFdF  dF ’  ~ijkl  dFijdFki  dFkl 

A  relationship  is  needed  between  stress  and  strain.  From  equation  (126),  one  can  see  that 


in  hyperelastic  materials,  P  is  related  to  F  through 

P  =  C  :  F,  4  =  CijkiFklt 

where 

F  =  du/dX  =  dv/dX , 
and  where  u  =  v  denotes  the  velocity. 

Substituting  equation  (127)  into  (124)  and  using  (128)  yields 

[  C  ::  (V05u  0  V0v)  dV=  f  t0  •  6udA  +  [  f0  •  6udV, 

Jv  JdV  Jv 

and  the  equivalent  indicial  form, 

f  f  M«,C IA+  f  faSmdV. 

jv  oAj  oJii  Jqy  Jv 


(127) 


(128) 


V<5u, 


(129) 


(130) 


This  is  the  virtual  work  equation  associated  with  hyperelasticity.  The  two-scale  approach  is 
described  next.  It  is  devised  so  that  traditional  finite  element  continuum  equations  can  be 
solved  in  the  coarse  scale  and  atomistic  equations  can  be  solved  in  the  fine  scale. 


3.3  Homogenization 


The  homogenization  framework  enables  the  weak  coupling  of  the  continuum  to  the  atoms. 
By  taking  the  limit  of  the  time-independent  asymptotic  expansion  parameter  e  -4  0,  we 
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exploit  the  weak  convergence  properties  of  the  scheme  so  as  to  decouple  the  length  scales. 
Hence  the  term  “weak  coupling.”  At  the  fine  scale,  the  domain  contains  only  atoms  with 
periodic  conditions  prescribed  on  the  boundary  and  all  atom  displacements  are  measured 
relative  to  a  fixed  point  in  the  local  frame  of  reference. 


The  homogenization  method  is  based  on  the  assumption  that  two  scales  exist  -  a  coarse 
scale  and  a  fine  scale.  Coordinates  in  the  coarse  material  scale  are  X  =  (Xi,X2,  A3),  and 
those  in  the  fine  material  scale  are  Y  =  (yi,y2,  Y3).  Likewise,  the  spatial  coordinates  are 
the  lowercase  analogues.  The  two  scales  are  related  by  the  scale  parameter 

Y  =  — .  (131) 

£ 

Therefore,  we  assume  that  the  ratio  of  scales  remains  the  same  before  and  after  deformation. 
The  aim  is  to  obtain  two  sets  of  coupled  equations.  The  asymptotic  series  assumption 
decomposes  the  displacements  as 


u(X)  =  u[0](X)  +  u[1](X) 


(132) 


=  u[0](X)+euW(Y), 


(133) 


where  ui°J  represents  the  displacement  at  the  coarse  scale  and  represents  the  perturbed 
displacements  due  to  inhomogeneity  at  the  fine  scale.  Square  brackets  denote  the  order  of 
the  term  in  the  asymptotic  series.  The  representation  of  the  total  displacement  at  the  fine 
scale  is  given  by  Takano  et  al.  [76]  as 


-u(X)  =  umicro(Y) 

€ 

=  F(u[°J(X))Y  +  uM(Y). 


(134) 


The  variable  X  in  equation  (134)  is  a  fixed  value  with  respect  to  Y.  That  is,  the  deformation 
gradient  of  a  point  in  the  coarse  scale  gets  mapped  onto  a  fine  scale  grid.  This  point  is 
typically  a  quadrature  point  in  a  finite  element  sense. 
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The  time  derivatives  are  analogous  to  equations  (133)  and  (134).  They  are  given  as 


u(X)  =  v(X) 

=  vI°1(X)  +  evW(Y), 

(135) 

micro  _  Vm^Cr°(X) 

=  F(v[0](X))Y  +  v[1](Y). 

(136) 

Substituting  equations  (133)  and  (135)  into  (130)  yields 

[  C  ::  [V*  (W°1(X)  +  <rfuW(Y))  0V*  (v™(X)  +evW(Y))]  dV 

Jv 

=  [  (<5ut0](X)+e<W1](Y))  -todA 

JdV 

+  f  (<5ut°](X)  +  etfuW(Y))  •  fo dV,  vW0],  <5uW 
Jv 


(137) 


Note  that  by  use  of  the  chain  rule  and  equation  (131), 

Vx<KX,Y)  =  V^+||vy0 

=  +  —  Vy0. 


(138) 


Therefore, 


Vx  (u[0>(X)  +  euW(Y))  =  Vxu«(X)  +  Vyu[1,(Y). 


(139) 


Using  equation  (139)  in  (137)  and  taking  the  average  over  Y  gives 

[  i  f  Cv.  [(Vx«uM(X)  +  VyiuW(Y))  ®  (VxvI°l(X)  +  VyvW(Y))]  dYdV 
Jv  \Y\  Jy 

=  f  («5U[°I(X)  +  £<5uM(Y))  •  todA  +  [  (W°](X)  +  eduW(Y))  •  f0dV,  (140) 

Jdv  Jv 

V<5u^, 

Then,  in  the  limit  as  e  -4  0,  equation  (140)  is  satisfied  only  if  the  following  two  equations 
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are  satisfied, 


if  [  C  ::  [Vx<5uM(X)  ®  (Vxvl°>(X)  +  VyvW(Y))]  dYdV 
Ml  Jv  Jy 

=  [  6uW(X)-i0dA+  [  <W°](X)  •  f0dV,  V  Su®, 

i/v’ 


(141) 


-L  f  [  C  ::  ;Vy<5ulll(Y)  0  (VxV^X)  -  VyvW(Y))]  rfl VJV  =  0,  VJu™.  (142) 

mi  Jv  Jy 

By  recourse  to  the  finite  element  method,  the  solution  of  equation  (141)  is  straightforward 
assuming  C  and  vM  are  known.  It  is  then  evident  that  due  to  the  dependence  of  (142)  on 
vt°l,  equations  (141)  and  (142)  are  coupled  and  must  be  solved  concurrently. 


In  the  next  two  sections,  a  method  is  first  shown  for  solving  equation  (142)  for  vM,  then  in 
the  following  section,  the  formulation  that  enables  the  atomistic  information  to  be  fed  into 
equation  (141)  is  derived.  Then  by  linearizing  the  equations,  a  Newton-Raphson  scheme  can 
be  employed  to  achieve  the  needed  concurrency. 


3.4  Atomistic  Equation 


Distinct  and  distinguishable  atoms  are  assumed  to  reside  in  the  local  level  cell.  By  the 
Cauchy-Born  rule  [12],  at  a  point  X,  F(u[oi)  is  assumed  to  give  the  energy  minimizing 
configuration  of  the  atoms.  For  simplicity,  we  assume  that  the  atoms  are  arranged  in  a 
lattice.*  Then,  the  positions  of  the  atoms  Y  are  given  from  the  lattice  coordinates  m  by 

Y(m)  =  me* :  m  6  £,£  =  Z3,Z  <  N,  (143) 

where  e*  are  the  primitive  translation  vectors  and  N  is  the  integer  multiple  of  atoms  con¬ 
tained  in  the  unit  cell.  To  avoid  confusion  in  notation,  atom  labels  are  noted  in  parentheses 

*Note  that  there  is  no  restriction  to  perfect  lattices.  In  fact,  by  use  of  computers,  arbitrary  arrangements 
of  atoms  can  be  considered  as  long  as  the  assumption  of  the  Cauchy-Born  rule  still  applies. 


99 


henceforth  and  are  not  subject  to  the  conventional  summation  rules  associated  with  indicial 
notation.  The  displacement  of  the  atoms  are 

q(m)  :  m  £  C.  (144) 

Upon  deformation,  the  new  positions  of  the  atoms  are  given  by 

y(m)  =  Y(m)  4-  q(m).  (145) 


The  deformation  gradient  is  defined  by 


F  = 


dy_ 

dY' 


(146) 


The  vector  separating  two  atoms  i  and  j  in  the  reference  configuration  is  given  by 


R(y)  —  Y(i)  Y(j), 


(147) 


where  Yy)  denotes  the  position  of  atom  j  and  Y(*)  the  position  of  atom  i.  The  vector 
separating  two  atoms  in  the  deformed  configuration  is  given  by 

rm  =  yu)-y  «•  (148) 


Then  the  Cauchy-Born  rule  can  be  stated  in  a  more  precise  manner  by 


T(ij)  =  FY0)  “  FY« 
=  FR(y). 


(149) 


For  the  energy  associated  with  the  deformation  of  the  atoms,  we  use  the  so-called  type  II 
parameterization  of  the  Tersoff-Brenner  potential  [94,95].  It  takes  the  form, 

W  =  -[£6(Y  +  q)-£6(Y)],  (150) 

na 

where  W  is  the  energy  density  of  the  frozen  system,  na  is  the  number  of  atoms,  and  E(,  is 
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the  binding  energy  given  for  a  pure  carbon  system  by 


=  E  E  lv*  (r<«>)  - Sv*  (r<«>)]  • 

*  i(>») 

&  =  2  -^yo)  > 

VpM  -  Mr)D(e)  -V5S/?(r-fi<«>) 

Rl  J  “  (5-1) 

K.w  = 

f  1,  rcflM 


/«)(»■)  = 


B{ij)  — 


I  1 1  +  cos 

0, 


r  >  J?<2> 


1  +  }  ]  G(9(ijk))f(ik)  ir(ik)) 


k(¥*J) 


—S 


G{e)  ~  °',{1  +  |“df 

with  the  constants  given  in  Table  3. 


+  (1  +  cos#)2 


J 

}' 


Table  3.  Parameters  for  Tersoff-Brenner  potential. 


1.39  A 

D (e) 

6.0  eV 

5 

1.22 

P 

2.1  A 

6 

0.5 

RW 

1.7  A 

bP> 

2.0  A 

a0 

0.00020813 

cl 

3302 

dl 

3.52 

(151) 

(152) 

(153) 

(154) 

(155) 

(156) 

(157) 


Given  that  the  energy  can  be  written  as  a  function  of  the  atom  displacements,  equation 
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(142)  can  be  expressed  in  a  form  conducive  to  atom  representations.  We  equate  vW  to  the 
rate  of  atom  displacement  and  attempt  to  solve  the  equivalent  form 


d  dvP  dCijki  dvf 
dYjijkl  dYi  ~  dYj  dXi' 

under  periodic  boundary  conditions.  The  solution  to  equation 
<97 Z  in  the  equation 


(158) 

(158)  is  found  as  the  zero  of 


am  =  /CvN  -  v  ■  v0vi°i(x), 


(159) 


where  K  is  the  na  x  na  Hessian  and  is  given  by 


K  = 


<92W 

dqdqj 


(160) 


where  q  is  the  vector  of  atom  displacements  of  size  3na  (in  three  dimensions),  and  is  a  third 
order  unsymmetric  tensor  that  is  obtained  from  the  first  derivative  of  the  Euler-Lagrange 
equation  with  respect  to  the  local  deformation  gradient,  given  by 


T>  = 


<92W 

dqdF' 


(161) 


The  size  of  depends  on  the  dimensionality  of  the  problem.  In  three  dimensions,  it  can  be 
expressed  as  an  na  x  9  matrix  where  9  corresponds  to  the  number  of  independent  components 
of  F. 


Under  the  Cauchy-Born  hypothesis  that  the  atomic  primitive  vectors  deform  homogeneously 
according  to  the  gross  deformation  of  the  continuum,  equation  (159)  is  a  linear  equation. 
That  is,  the  atom  positions  needed  to  compute  K  are  determined  solely  by  Vov^  and  is 
independent  of  v^.*  However,  if  necessary,  a  mixed  strategy  can  be  used  incorporating 
molecular  dynamics  to  compute  the  true  energy  minimizing  atom  configuration  at  the  fine 
scale. 

•This  is  not  to  say  that  v^l  and  v—  are  independent.  The  distinction  being  made  here  is  between 
interscale  and  intrascale  coupling.  The  coupling  is  interscale  but  not  intrascale. 
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3.5  Multiscale  Equation 


Once  equation  (159)  has  been  solved  for  the  remaining  task  is  to  formulate  a  tractable 
global  scale  boundary  value  problem.  The  key  distinction  between  this  investigation  and 
conventional  continuum  formulations,  such  as  hyperelasticity,  is  the  conspicuous  incorpora¬ 
tion  of  vM,  a  fine-scale/atomistic  quantity,  in  the  global  scale  equations,  and  the  definition 
of  the  material  property  tensor  completely  in  terms  of  atomistic  variables. 


We  return  to  equation  (141)  recognizing  that  vW  is  now  known.  Incorporating  the  definition 
for  the  first  Lagrangian  elasticity  tensor  from  equation  (126)  yields 


u 


\Y\JvJy  dFdF 


[V*<W0](X)  <g>  (VXV[0](X)  +  VyvW(Y))]  dYdV 

[  6u[0]{X)-t0dA+  [  <W°](X)-f0dV,  V  Sn[0]. 
JdV  Jv 


(162) 


Then  using  the  definition  of  F  in  equation  (118)  and  assuming  a  first  order  Taylor  series 
representation  for  the  time  derivative  gives, 


d2W 

dFdF 


::  (Vx<5uM(X)  ®  V*V[01(X))  dYdV 


=  [  <Su[01(X)  •  to dA  +  [  du[0](X)  •  f0dV 
JdV  Jv 


dFdq 


::  (Vxtfu[°l(X)  0  vW(Y))  dYdV,  V  <5u[o1. 


The  solution  to  equation  (163)  yields  that  takes  into  consideration  the  effect  of  the  atoms. 
It  is  noteworthy  that  the  last  term  is  zero  when  the  energy  distribution  over  Y  is  constant, 
i.e.,  when  the  atom  arrangement  forms  a  perfect  lattice.  This  reduces  the  problem  to  a 
classical  harmonic  approximation  where  the  first  Lagrangian  elasticity  tensor  is  assumed  to 
model  the  material  behavior.  In  equation  (163),  the  last  term  serves  as  a  corrective  force  in 
regions  of  highly  energetic  atoms,  i.e.  non-locality  regions,  to  account  for  defects  and  lattice 
inhomogeneities.  Note  also  that  the  atomistic  energy  density  and  its  derivatives  intrinsically 
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account  for  invariance  properties  of  the  atoms.  For  example,  by  switching  the  positions  of 
two  atoms  of  the  same  species,  the  energy  remains  constant. 


3.6  The  Euler- Lagrange  Equations  and  the  Hessian 


In  this  section,  the  analytic  forms  of  the  Euler-Lagrange  equations  and  Hessian  are  derived 
for  a  general  potential.  The  nontrivial  algebra  typically  needed  to  obtain  equations  (160) 
and  (161)  for  the  specific  case  of  the  Tersoff-Brenner  potential  are  shown  in  greater  detail  in 
the  Appendix,  and  only  general  forms  are  derived  here.  The  Euler-Lagrange  equation  is  the 
first  derivative  of  the  Lagrangian  with  respect  to  the  degrees  of  freedom.  In  this  problem, 
the  Lagrangian  is  the  negative  of  the  atomistic  energy  density.  The  Euler-Lagrange  equation 
is  therefore  given  by 


s  =  J>w_ 

dq  (m) 

=  1  dEb 
na  dq[m)  ’ 

and  using  the  chain  rule  for  derivatives,  it  is 

£=1  m 

na  ^q(ro) 

_  J_  /  dEb  dr{l]}  |  dEt  dr{lk)  |  dEb  drUk)  \ 

na  (9q(m)  dv^\  3q(m)  dr^j^  dq^m)  J 


(164) 


(165) 


The  Hessian  is  obtained  by  taking  an  additional  derivative  of  the  Euler-Lagrange  equations. 
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Specifically,  we  again  make  use  of  the  chain  rule  to  obtain 
d2W 


K 


dq(n)dq(m) 

1  &Eb 
Tla  ^q(n)^q(m) 
d2Eh 


n„ 


+ 


dr  ms  dr. 


..  ( ?Em  0  +  92 Eb 


{ij)OT{ij) 


d2Eb 


\dq  (n) 

(  drm 


<8> 


+ 


+ 


+ 


dq  (m) ) 
^  dT(lj)  ) 

dr(ik)dr{ij)  "  V  dq(n)  ^  dq (TO)/ 
d^Eb  ..  / 0  g£W 
drUk)dr{jk)  ^q(m) 

c>2£<,  ..  ( 


dT(jk)  0  ^r(jj)  ^ 


5r0fc)ar(y)  "  V  dq(n)  dq{m)  J 
d2Eb  „  ( dtm  Q  df(ifc)  Y 


+ 


^r(y)^r(ifc) 

a2^6 

d^[ih)d^[ik) 

d2Eb 

d*(ij)dT(jk) 

d2Eb 

dr(ih)dr(jk) 


dr{jk)dv{ik)  "  V^q(n)  5q(m)/J  ’ 

Second  derivatives  of  the  inter-atom  vectors  are  zero,  i.e., 
d2r{ij)  _  d2 


dr  ft)  dr 


<8> 


(**) 


"  V^q(n)  9q(m)> 
/  dr(ik)  Q  gr^)  \ 

V^q(n)  dq(m)  / 

.  f  arfa)  0  \ 

'  V^q(n)  dqim)J 

( dr(ik)  0  drw 
V9q(n)  dq(rn) 


dq{m)dq{n) 


dq{m)dq(n) 

=  0  V(m,  n). 


(Y(j)  +  <Jo>  -  Y(j)  +  qiii) 


(166) 


(167) 

(168) 


Next,  the  appropriate  right  hand  side  expressions  are  derived  for  equation  (159)  and  (161). 
This  involves  the  use  of  the  chain  rule  again  to  obtain 


d2W  _  2 

dq{m)dF  na 


d2Eb 


\.d*(ij)dr(tj) 

d2Eb 

dr^ik'jdr^ij'j 

d2Eb 


dr<fj)drm 

d2Eb 


and  by  definition, 


dr 


(*J) 


d  F 

dT{ij) 

dF 

drm 


dF 

_  /  dr<jk) 

dr{ik)drijk)  "  \  dF 


i+ 

^Q(m)  J 

\  + 

1 + 

®  ar<«) ' 
^l(m)  > 

)+ 

aq(m)  y 

1+ 

d2Eb 


dr  (ij'jdr  (ik'j 

d2Eb 

dv(ik)dr{ik ) 
d2Eb 

dr{3k)dr{lj) 

d2Eb 

dr{jk)dr[ik) 

d2Eb 


dr  (i 


( ik ) 


dr, 


Xv) 


drUk)drUk)  "  \  dF  dq( 


( 


dF  dq  (TO) 
drm  dr{ik) 
dF  dq(m) 

aF  ^aq(m) 

( drm  drUk) 
\  dF  dq(m) 
drm  0  5r0,} 


^r(d)  _  -p  dr(jk)  _ 


dF 


dF 


=  Rfi 


'(**)• 


(169) 


(170) 
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Finally,  we  use  a  similar  approach  to  define  the  first  Lagrangian  elasticity  tensor.  This  is 
the  traditional  way  of  estimating  the  elastic  properties  of  a  solid.  Using  the  chain  rule  once 
again  gives 


d2W 

dFdF 


1_ 


d2Eb 


dr{ij)dr{ij) 
dT{ik)dr(ij) 

]  &Eb 
dr(ij)dT(jk) 

| 

^{ik)^(jk) 


d2Eb 

dr{ij)dT{ik) 

d2Eb 

dv{ik)dr{ik) 

d2Eb 

dr{jk)dr{ij) 

d2Eb 

dr(jk)dr(ik) 
d2Eb 


(171) 


d*Uk)drUk)  "  \  dF  v  dF 

The  first  Lagrangian  elasticity  tensor  is  used  in  equation  (163)  whose  solution  gives  vt°l. 


In 


a  perfect  lattice,  equation  (171)  provides  the  only  atomistic  material  information  needed  to 
solve  the  macroscopic  continuum  problem.  The  next  section  illustrates  this  by  showing  that 
the  perturbation  is  zero  for  a  uniform  crystal. 


3.7  Example  Problems 

3.7.1  Example  I:  Perfect  1-D  Atomic  Lattice 

To  illustrate  the  calculation,  a  1-D  analytical  example  is  presented.  The  TersofF-Brenner 
potential  is  used  to  represent  the  energetics  of  a  1-D  single-species  chain  of  carbon  atoms. 
The  objective  here  is  to  solve  equations  (158)  and  (159)  for  and  demonstrate  a  simple 
case  of  a  perfect  lattice  using  this  method. 

One  atom  comprises  the  periodic  unit  cell,  but  to  account  for  the  effects  of  triples,  two 
“fictitious”  atoms  are  assumed  to  extend  beyond  the  boundaries  of  the  cell  on  each  side  as 
illustrated  in  Figure  28.  Periodic  conditions  apply  at  the  cell  boundaries.  The  equilibrium 
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Fictitious  atoms 

cf  \) 


Fictitious  atoms 

/  \ 

o  o 


\  / 

Cell  boundaries 

Figure  28.  Unit  cell  of  1-D  carbon  chain.  The  atoms  are  labeled  by  identify¬ 
ing  numbers. 


lattice  constant  for  the  one-dimensional  chain  is  r0  =  1.86868A  The  following  two  conditions 
stem  from  the  1-D  assumption, 


8  =  7T, 

Rw  <  r  <  i?(2). 


(172) 

(173) 


This  simplifies  the  expressions  in  the  Appendix.  The  resulting  Hessian  for  three  arbitrary 
collinear  atoms  (i.  j,  k)  is  obtained  as 


r-m) 

*'11 

jy-iijk) 

'''12 

y{ijk) 

*'13 

[JC{ijk)]  = 

^(yfc) 

Ml 

yiijk) 

”'22 

r(ijfc) 

”'23 

y-iijk) 

''■'31 

y(ijk) 

^32 

r(bfc) 

*'33 

(174) 


where  =  Knm' )  and  the  terms  are  defined  by 

r<f  >  =VR  -  BV :  +  a06V'AB'+*  f[ik)  +  \(S  + 1) VAB$  (ajmy 

.  a°^\ /  D1  +  J  f" 


«Sf 1 = -  v* + bv;  -  4>. 

f  >  =  _  ^V-AB$/m  -  5 (*  +  1)VaB$  -  a-£vA B$ ft 

K%k)  =v£  -  BVL 


jr(b 

^13 


'ifc)  5 


-22  “ 


r{ijk)  Da+7  f' 

^23  “  /(ifc)j 


4f ’  =4(<S  +  1)^,1?$  (<•„/('«, )2  +  ^VAB$fm- 


(175) 

(176) 

(177) 

(178) 

(179) 

(180) 
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Upon  assembly  of  the  two  unique  pairs  (31,12)  and  their  associated  triples  (123,132)  in 
Figure  28,  the  final  assembled  Hessian  of  the  global  system  is  given  by  the  matrix, 


[JC}  = 


JCu 

^13 

X-21 

^22 

£23 

^31 

£32 

K33 

which  is  assembled  through  the  operation, 


(181) 


[/C]=[_JU  =  /Cm„,  (182) 

(m)  (n) 

where  |J  is  the  addition  operator  over  all  unique  pair  and  triple  combinations  of  (i,j,  k)  and 
(m)  and  (n)  are  displacement  degrees  of  freedom  for  each  atom.  In  equation  (182),  [/C]  is 
symmetric  once  again  and  its  components  are  obtained  in  detail  for  the  problem  shown  in 
Figure  28  as  follows, 


Kn  =  4123) + 43215)  +  4?41) ,  (183) 

=  /c£23)  +  /Cg41),  (184) 

/C13  =  /cS323)+/cS3215),  (185) 

fC22  =  41223)  +  ;cn41).  (186) 

/c23  =  /cg23),  (187) 

JC33  =  /Cg23)  +  /CS315).  (188) 


This  constitutes  the  stiffness  matrix  /C  in  equation  (159). 


The  next  step  is  to  calculate  the  right-hand  side  of  equation  (158)  which  is  equivalent  to  cal¬ 
culating  T>  and  multiplying  by  the  global  rate  of  the  deformation  gradient.  Using  equations 
(161)  and  (169)  the  right  hand  side  for  three  arbitrary  collinear  atoms  (i,j,k)  is  obtained  as 

'  v[ijk) 

<  v{ijk) 

'ni.ijk) 

l  ^3 
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where  the  components  are  defined  by 


©p  =Rlm  (v;  -  BV")  +  {Riii)  -  Rm)  (^Xfl  A>) 

+  Rm  (<s  + 1 )vaB{v)"  (®<>/(it))  +  ~YvAB(ij‘ f«k)j  > 

=  -  %i)  (K  -  BVl)  -  Rm  (^vX«  4>) , 

=  -  ■%>  (^vX<5>%))  ~  Rm  (j(«  +  VVaB$  (oXi)’  + 


(190) 


(192) 


As  earlier,  the  assembly  operation 

{V}=  □  {^}=Pm 

(m) 

yields  the  right-hand  side  of  the  global  system  given  by, 


where  the  components  are 


(193) 


(194) 


K1=P<I23)+pf5>+P<24I),  (195) 

X>2  =  b'123)  +D;241),  (196) 

Z>,  =  V'3la>  +  ®<315).  (197) 


Under  the  assumption  of  a  1-D  perfect  lattice,  we  have  =  R(ik),  and  consequently, 
V\  =  0.  Then,  we  can  satisfy  the  periodicity  condition  and  the  rigid  body  constraint  by 
setting  njgj  =  =  0.  The  solution  is  therefore, 


-  «W  - 

;(i) 


=  v , 


(2)  ~  V(3) 


0. 


(198) 


In  light  of  equation  (198),  the  last  term  in  equation  (163)  is  zero  and  the  material  properties 
are  obtained  from  the  atomistic  energy  density  solely  through  equation  (171).  This  result 
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shows  that  in  a  defect-free  lattice,  the  homogenization  method  coincides  with  the  conven¬ 
tional  atomistic  hyperelasticity  problem.  The  next  section  shows  an  example  in  which  a 
defect  causes  an  inhomogeneous  energy  distribution  leading  to  a  situation  where  homoge¬ 
nization  is  needed  to  average  out  the  energy. 


3.7.2  Example  II:  1-D  Atomic  Lattice  With  Defect 


Consider  the  problem  shown  in  Figure  29  where  the  center  atom  is  displaced  by  a  distance 
L  from  its  original  energy  minimizing  configuration.  This  displacement  of  the  center  atom 
constitutes  the  defect.  With  this  change,  the  assumption  in  equation  (173)  no  longer  applies 
and  the  key  stiffness  matrix  term  in  equation  (183)  is  now 


Ku  ~VRn-n  ~  BWVA[12)  +  a°SVAil2)B(12)  /(13) 


*(12) 

s 


+  5<*  +  i  )va^b\U 


(g°/(13))  +  ~YVA{l2)B\u)  f (13) 


+  Vr(1S)  ”  B(13)VA(1 3)  +  a°SVAil3)B(13)  / (12) 

+  \  (*  +  1  .)W$  («./,', 2,)2  +  ^  K‘+ J 


(199) 


2  '  ^(13)  B(13)f(12)  > 

and  likewise,  equation  (195)  becomes, 

T>i  =  72(12)  (Vr(12)  -  B{i2)VA (1J))  +  (72(13)  -  72(i2))  (jyva{12)B(i2)  /(i3) 

+  72(13)  ^2^  +  1)Va(12)£(i2)  (ao/(13))  +“|"Va(i2)-B(12)4/(X3)) 

-  #(13)  (Vr(13)  “  #(13)  ^(13))  -  (#(12)  -  #(13))  (^FAa3)#(\tf/(12)) 

-  #(12)  +  1  )VA(13)Blu)  (°o/(  12))  +  ^^(13)  #(13?  /(12))  • 

where  72(i2)  =  r0  -  L  and  72(i3)  =  r0  +  L.  Then,  solving  equation  (159)  under  periodic 
boundary  conditions  gives 


(200) 


uf]  =  ^V0v(°], 


=  yW  _  o. 


(201) 
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Figure  29.  Unit  cell  of  1-D  carbon  chain  with  periodic  defect. 


The  /Vo^°l  solution  as  a  function  of  L/r0  is  shown  in  Figure  30.  As  expected,  the  solution 
has  symmetry  about  the  origin  and  grows  asymptotically  larger  as  the  size  of  the  defect  ( L ) 
grows  closer  to  the  cut-off  radii.  Larger  defects  are  avoided  presently  due  to  the  non-convex 
structure  of  the  energy  well  associated  with  the  Tersoff-Brenner  potential.  This  generally 
leads  to  unphysical  discontinuities  in  the  perturbation  velocity  (v^)  due  to  discontinuous 
second  derivatives  of  the  atomistic  energy  with  respect  to  the  defect  size.  This  is  attributable 
to  the  construction  of  the  empirical  potential  in  equations  (151)— (157)  which  is  intrinsically 
suited  for  systems  where  nearest  neighbor  atoms,  even  in  defect  regions,  are  within  the  cut-off 
radius  Rf®. 

It  is  also  noteworthy  that  arbitrary  defect  densities  can  be  treated  by  appropriate  modifica¬ 
tion  of  the  unit  cell.  In  most  cases,  one  can  tailor  the  desired  density  by  increasing  the  size 
of  the  unit  cell  and  performing  the  summations  and  the  assembly  of  the  atomistic  discrete 
equations  over  more  atoms.  Figure  31  illustrates  this  idea  for  the  1-D  carbon  chain. 

Numerical  experiments  show  that  as  the  size  of  the  unit  cell  increases,  the  perturbative 
displacement  has  a  sharp  discontinuity  at  the  defect.  Figure  32  shows  this  non-local  behavior 
as  the  number  of  atoms  increases.  The  problem  is  of  a  single  defect  in  chains  of  increasing 
size.  The  defect  magnitude  is  held  fixed  at  L/r0  =  0.01.  The  non-local  discontinuity  of  the 
perturbative  velocity  qualitatively  agrees  with  traditional  displacement  jumps  that  occur  at 
dislocation  cores.  The  discontinuity  indicates  that  the  material  property  at  the  defect  ( Jp||0 
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Figure  30.  Distribution  of  t>M /V0u^  solution  as  a  function  of  the  defect  size. 

O  O  i  ©  ©  ©  ©  ©  •©©©©©:  O  O 

;  \ 

Defect 

Figure  31.  Larger  chain  of  atoms  in  perfect  arrangement  around  the  defect 
region  decreases  the  defect  density. 

is  modified  by  the  last  term  in  equation  (163),  an  amount  proportional  to  u[1]  that  serves  as 
a  correcting  force  for  the  non-locality. 

Although  the  primary  details  of  the  method  have  been  demonstrated  in  these  two  examples, 
the  method  can  be  extended  to  consider  the  multiscale  problem  shown  in  equation  (163) 
for  more  general  cases  involving  self-consistent  solutions  with  equation  (159).  We  presently 
restrict  ourselves  to  analytical  1-D  examples  where  the  results,  as  in  the  previous  examples, 
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Figure  32.  Distribution  of  uW/V0u^  along  unit  cell  length  for  varying  num¬ 
ber  of  atoms  ( L/r0  =  0.01). 

can  be  reported  independent  of  the  macroscale  solution  v^°\  and  reserve  more  complex  cases 
in  higher  dimensions  for  a  separate  numerical  investigation. 


3.8  Closing  Remarks 


Linking  atomic  scale  physics  with  continuum  scale  phenomena  is  of  keen  interest  in  the  study 
of  failure,  fracture,  and  reliability  of  engineering  structures.  The  effects  that  dominate  the 
failure  at  the  continuum  scale  typically  initiate  and  evolve  from  the  atomic  scale.  Despite 
numerous  promising  methods  in  the  literature  that  are  capable  of  linking  scales  up  to  the 
micron  level,  structures  at  the  meter  scale  and  beyond  have  only  begun  to  be  studied.  To  this 
end,  in  this  part  we  have  attempted  to  address  this  issue  by  exploiting  the  weak  convergence 
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properties  of  homogenization  to  devise  a  scheme  which  passes  atomistic  information  to  very- 
large  continuum  scales. 

We  have  applied  the  Cauchy  Born  rule  [12]  to  the  atom  scale  by  assuming  that  the  config¬ 
uration  of  atoms  used  to  solve  for  the  perturbation  displacement  is  indeed  the  minimizing 
configuration  of  the  atomistic  energy.  We  have  not  considered  the  method  in  conjunction 
with  a  molecular  dynamics  routine,  i.e.,  various  strategies  of  minimizing  the  atomistic  en¬ 
ergy  by  quenching  through  artificial  temperature  decrease  or  solving  Newton’s  equations  to 
minimize  the  interatom  forces. 

For  this  work,  the  specific  case  of  the  Tersoff-Brenner  type  II  potential  was  considered. 
But  the  principles  and  the  general  equations  can  be  extended  to  any  potential  provided  the 
appropriate  derivatives  can  be  obtained  as  in  Appendix  A.  Typically  for  classical  systems, 
onerous  tensor  algebra  and  calculus  are  required. 

The  aim  of  this  part  was  to  develop  an  approach  by  which  atomistic  physics  can  be  em¬ 
bedded  into  a  continuum  formulation  for  large  scale  systems.  This  goal  has  been  achieved 
by  formulating  a  consistent  set  of  equations  involving  a  classical  atomistic  potential  at  the 
fine  scale  and  general  finite  strain  and  deformation  elasticity  at  the  coarse  scale.  Simple  1-D 
analytical  results  were  shown  to  illustrate  the  approach  and  its  features.  More  realistic  mul- 
tiaxial  problems  in  two  and  three  dimensions  for  more  detailed  validation  are  the  subjects 
of  ongoing  work. 
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Appendix  A.  Derivatives  of  the 
Tersoff- Brenner  Potential 


The  derivatives  needed  to  form  the  Euler-Lagrange  equations  and  the  Hessian  are  shown 
here  in  detail.  To  simplify  the  notation,  we  define  the  following  expressions, 

r(y)  =  |r(y)|  >  (A-l) 
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Note  that  although  the  equations  are  written  in  component  form  with  respect  to  atoms,  it 
is  still  in  dyadic  notation  due  to  the  multi-axial  components  of  That  is  r^)  •  ex  is  the 
component  of  the  vector  originating  at  atom  i  and  terminating  at  atom  j  in  the  direction 
of  ex,  r(jj)  •  e2  is  the  component  in  the  direction  of  e2,  etc.  Referring  to  equations  in  this 
report,  (151)— (157) ,  the  derivatives  in  equation  (165)  are  defined  by 
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for  (mn)  =  (z'j),  (ife)  when  7  =  (ijfc)  and  (ran)  =  (ij),  (;'&)  when  7  =  (jz'fc).  The  angles  %*) 

and  %zfc),  shown  in  Figure  A-l,  are  the  angles  subtending  the  connecting  lines  at  the  atoms 

i  and  j,  respectively.  Note  that  r^)  =  - r y<).  The  following  identities  can  also  be  shown: 
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Figure  A-l.  Angles  and  interatom  vectors. 


where  I  is  the  3x3  identity  tensor  for  a  system  in  a  3-D  domain.  The  derivative  with  respect 
to  r(jfc)  can  be  obtained  likewise.  This  indicates  that  the  summations  in  equations  (A-3), 
(A-4),  (A-8),  and  (A-9),  when  multiplied  by  equation  (A-17)  in  equation  (165)  are  nontrivial 
if  and  only  if  m  is  equal  to  i,  j,  or  k. 


For  the  Hessian  in  equation  (166),  the  differential  terms  are  defined  by, 
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To  complete  the  derivation,  the  following  identities  are  needed, 
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and  ( mn,pq )  =  ( ij,jk ),  ( jk,ij ),  ( jkjk )  for  7  =  jik,  and, 
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This  completes  the  closed-form  derivation  of  the  Hessian  for  the  Tersoff-Brenner1,2  potential. 


Despite  the  relative  algebraic  complexity  of  the  expressions,  the  calculations  can  be  per¬ 
formed  readily  using  computers.  The  algorithm  is  based  on  an  additive  assembly  process 
by  casting  the  equations  in  their  equivalent  matrix  forms  and  then  summing  over  all  unique 
pairs  and  triples  of  atoms,  which  translates  well  to  an  iterative  computational  methodology. 


1Tersoff,  J.  “Empirical  Interatomic  Potential  for  Carbon,  With  Applications  to  Amorphous  Carbon.57 
Physical  Review  Letters ,  voh  61,  no.  25,  pp.  2879-2882,  19  December  1988. 

2 Brenner,  D.  W.  “Empirical  Potential  for  Hydrocarbons  for  Use  in  Simulating  Chemical  Vapor  Deposition 
of  Diamond  Films.55  Physical  Review  B,  vol.  42,  no.  15,  pp.  9458-9471,  15  November  1990. 
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